{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "4dbd5617-33f8-4544-82d3-0f96e82c50e6",
    "showInput": false
   },
   "source": [
    "## Tutorial: **Hidden Markov model**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "f699d5cc-9798-4828-8732-70d2d497edba",
    "showInput": false
   },
   "source": [
    "This tutorial demonstrates modeling and running inference on a hidden Markov model (HMM)\n",
    "in Bean Machine. The flexibility of this model allows us to demonstrate some of the\n",
    "great unique features of Bean Machine, such as block inference, compositional inference,\n",
    "and separation of data from the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "9390c28c-9c48-48cc-a207-3f73dbabbcc6",
    "showInput": false,
    "tags": []
   },
   "source": [
    "## Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "9390c28c-9c48-48cc-a207-3f73dbabbcc6",
    "showInput": false
   },
   "source": [
    "HMMs are a class of probabilistic models which are popular for doing inference on\n",
    "discrete-time stochastic processes. In general, Markov models are used to study a\n",
    "sequence of random variables, $X_1,\\ldots,X_N$, where the sequence is \"memoryless\" such\n",
    "that the distribution of $X_{n}$ depends only on the value of $X_{n-1}$; any sequence\n",
    "which is memoryless is said to satisfy the Markov property. One reason Markov models are\n",
    "popular is because this flexible framework can be used to model many time-evolving\n",
    "processes, such as words in a sentence, the position of a robot, or the weather.\n",
    "\n",
    "An HMM is a Markov model in which observations are modeled as being noisy. While we are\n",
    "interested in doing inference on each $X_n$, we are actually observing variables\n",
    "$Y_1,\\ldots,Y_N$ which can depend on the values of $X_1,\\ldots,X_N$ in a variety of\n",
    "ways. In specific settings, HMMs can be very tractable, and lend themselves towards\n",
    "provably-optimal algorithms such as Kalman Filtering. Here, we illustrate how to do\n",
    "general inference with MCMC as applicable to general HMMs. The single-site algorithms\n",
    "underlying Bean Machine enable inference algorithms which scale favorably with the size\n",
    "of the HMM."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "a50a9e55-266a-4e70-be95-e7cc18e8381e",
    "showInput": false
   },
   "source": [
    "## Prerequisites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "d5858da4-abf4-47ad-9cf3-f2c415faeb0e",
    "showInput": false
   },
   "source": [
    "Let's code this in Bean Machine! Import the Bean Machine library, some fundamental PyTorch classes, and optionally typing for our code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install Bean Machine in Colab if using Colab.\n",
    "import sys\n",
    "\n",
    "\n",
    "if \"google.colab\" in sys.modules and \"beanmachine\" not in sys.modules:\n",
    "    !pip install beanmachine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import os\n",
    "import warnings\n",
    "\n",
    "import beanmachine.ppl as bm\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "import torch.distributions as dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell includes convenient configuration settings to improve the notebook\n",
    "presentation as well as setting a manual seed for reproducibility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Eliminate excess inference logging from Bean Machine, except for progress bars.\n",
    "logging.getLogger(\"beanmachine\").setLevel(50)\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Manual seed\n",
    "bm.seed(111)\n",
    "\n",
    "# Other settings for the notebook.\n",
    "smoke_test = \"SANDCASTLE_NEXUS\" in os.environ or \"CI\" in os.environ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "d340d033-17f1-423c-9734-823e3ee18062",
    "showInput": false
   },
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "d340d033-17f1-423c-9734-823e3ee18062",
    "showInput": false
   },
   "source": [
    "We model the hidden state $X$ as being a discrete-time Markov chain with $K$ states and\n",
    "transition matrix (_a.k.a._ kernel) $\\Theta$. We model $N$ time steps of this chain with\n",
    "variables $X_1,\\ldots,X_N$, and use variables $Y_1,\\ldots,Y_N$ to model observable\n",
    "emissions of each $X$ with Gaussian noise.\n",
    "\n",
    "Formally, the transition and emission probabilities are as follows, for\n",
    "$n\\in1,\\ldots,N$:\n",
    "\n",
    "- $X_{n+1}\\mid X_n\\sim\\text{Categorical}(\\Theta[X_n])$\n",
    "- $Y_n\\mid X_n\\sim\\text{Normal}(\\mu[X_n],\\sigma[X_n])$\n",
    "\n",
    "\n",
    "Accordingly, priors can be assigned as follows, for $k\\in1,\\ldots,K$:\n",
    "\n",
    "- $\\mu[k]\\sim\\text{Normal}(\\mu_\\text{loc},\\mu_\\text{scale})$\n",
    "- $\\sigma[k]\\sim\\text{Gamma}(\\sigma_\\text{shape},\\sigma_\\text{rate})$\n",
    "- $\\Theta[k]\\sim\\text{Dirichlet}([\\frac{c}{K},\\ldots,\\frac{c}{K}])$\n",
    "\n",
    "\n",
    "Finally, assume that the value of $X_1$ is known:\n",
    "\n",
    "- $X_1=0$\n",
    "\n",
    "\n",
    "So the model is set by choosing the prior through the values of:\n",
    "$\\mu_\\text{loc},\\mu_\\text{scale},\\sigma_\\text{shape},\\sigma_\\text{rate},c$. ($c$ stands\n",
    "for concentration).\n",
    "\n",
    "We can implement this model in Bean Machine by defining random variable objects with the\n",
    "`@bm.random_variable` decorator. These functions behave differently than ordinary Python\n",
    "functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "d340d033-17f1-423c-9734-823e3ee18062",
    "showInput": false
   },
   "source": [
    "<div\n",
    "  style=\"background: #daeaf3;\n",
    "      border_left: 3px solid #2980b9;\n",
    "      display: block;\n",
    "      margin: 16px 0;\n",
    "      padding: 12px;\"\n",
    ">\n",
    "  Semantics for <code>@bm.random_variable</code> functions:\n",
    "  <ul>\n",
    "    <li>\n",
    "      They must return PyTorch <code>Distribution</code> objects.\n",
    "    </li>\n",
    "    <li>\n",
    "      Though they return distributions, callees actually receive <i>samples</i> from the\n",
    "      distribution. The machinery for obtaining samples from distributions is handled\n",
    "      internally by Bean Machine.\n",
    "    </li>\n",
    "    <li>\n",
    "      Inference runs the model through many iterations. During a particular inference\n",
    "      iteration, a distinct random variable will correspond to exactly one sampled\n",
    "      value: <b>calls to the same random variable function with the same arguments will\n",
    "      receive the same sampled value within one inference iteration</b>. This makes it\n",
    "      easy for multiple components of your model to refer to the same logical random\n",
    "      variable.\n",
    "    </li>\n",
    "    <li>\n",
    "      Consequently, to define distinct random variables that correspond to different\n",
    "      sampled values during a particular inference iteration, an effective practice is\n",
    "      to add a dummy \"indexing\" parameter to the function. Distinct random variables\n",
    "      can be referred to with different values for this index.\n",
    "    </li>\n",
    "    <li>\n",
    "      Please see the documentation for more information about this decorator.\n",
    "    </li>\n",
    "  </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "d340d033-17f1-423c-9734-823e3ee18062",
    "showInput": false
   },
   "source": [
    "Note also that, compared to the statistical notation above, our implementation uses\n",
    "0-indexing instead of 1-indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HiddenMarkovModel:\n",
    "    def __init__(\n",
    "        self,\n",
    "        N: int,\n",
    "        K: int,\n",
    "        concentration: float,\n",
    "        mu_loc: float,\n",
    "        mu_scale: float,\n",
    "        sigma_shape: float,\n",
    "        sigma_rate: float,\n",
    "    ) -> None:\n",
    "        self.N = N\n",
    "        self.K = K\n",
    "\n",
    "        self.concentration = concentration\n",
    "        self.mu_loc = mu_loc\n",
    "        self.mu_scale = mu_scale\n",
    "        self.sigma_shape = sigma_shape\n",
    "        self.sigma_rate = sigma_rate\n",
    "\n",
    "    @bm.random_variable\n",
    "    def Theta(self, k):\n",
    "        return dist.Dirichlet(torch.ones(self.K) * self.concentration / self.K)\n",
    "\n",
    "    @bm.random_variable\n",
    "    def Mu(self, k):\n",
    "        return dist.Normal(self.mu_loc, self.mu_scale)\n",
    "\n",
    "    @bm.random_variable\n",
    "    def Sigma(self, k):\n",
    "        return dist.Gamma(self.sigma_shape, self.sigma_rate)\n",
    "\n",
    "    @bm.random_variable\n",
    "    def X(self, n: int):\n",
    "        if n == 0:\n",
    "            return dist.Categorical(torch.tensor([1.0] + [0.0] * (self.K - 1)))\n",
    "        else:\n",
    "            return dist.Categorical(self.Theta(self.X(n - 1).item()))\n",
    "\n",
    "    @bm.random_variable\n",
    "    def Y(self, n: int):\n",
    "        return dist.Normal(self.Mu(self.X(n).item()), self.Sigma(self.X(n).item()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will generate random observations and choose the priors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_chain_observations(theta, mus, sigmas, N):\n",
    "    theta_distbns = {j: dist.Categorical(vector) for j, vector in enumerate(theta)}\n",
    "    hidden = [0]\n",
    "    while len(hidden) < N:\n",
    "        hidden.append(theta_distbns[hidden[-1]].sample().item())\n",
    "\n",
    "    def observe(k):\n",
    "        return dist.Normal(mus[k], sigmas[k]).sample().item()\n",
    "\n",
    "    return hidden, list(map(observe, hidden))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "concentration = 1.1\n",
    "mu_loc = 0.0\n",
    "mu_scale = 5.0\n",
    "sigma_shape = 0.5\n",
    "sigma_rate = 1.0\n",
    "\n",
    "N = 15\n",
    "K = 2\n",
    "\n",
    "thetas_obs = dist.Dirichlet(torch.ones(K) * concentration / K).sample((K,))\n",
    "mus_obs = dist.Normal(mu_loc, mu_scale).sample((K,))\n",
    "sigmas_obs = dist.Gamma(sigma_shape, sigma_rate).sample((K,))\n",
    "\n",
    "\n",
    "xs_obs, ys_obs = generate_chain_observations(thetas_obs, mus_obs, sigmas_obs, N)\n",
    "x_obs = torch.tensor(xs_obs)\n",
    "y_obs = torch.tensor(ys_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize model\n",
    "model = HiddenMarkovModel(\n",
    "    N,\n",
    "    K,\n",
    "    concentration,\n",
    "    mu_loc,\n",
    "    mu_scale,\n",
    "    sigma_shape,\n",
    "    sigma_rate,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inference is the process of combining _model_ with _data_ to obtain _insights_, in the\n",
    "form of probability distributions over values of interest. Bean Machine offers a\n",
    "powerful and general inference framework to enable fitting arbitrary models to data.\n",
    "\n",
    "Our model makes use of both continuous and discrete random variables. We'll want to make\n",
    "use of different inference strategies for each. In particular, we would like to take\n",
    "advantage of gradient information for the continuous random variables. To facilitate\n",
    "this, Bean Machine provides the `CompositionalInference` class.\n",
    "\n",
    "`CompositionalInference` is a powerful, flexible class for configuring inference in a\n",
    "variety of ways. By default, `CompositionalInference` will select an inference method\n",
    "for each random variable that is appropriate based on its support. The HMM presented in\n",
    "this tutorial has a number of different random variables that we're interested in\n",
    "learning about, and we will customize our own. The proposers for each family of random variables are summarized in the table below:\n",
    "\n",
    "| Random variable | Support            | Inference method                            |\n",
    "| --------------- | ------------------ | ------------------------------------------- |\n",
    "| $X$             | $0,\\ldots,K-1$     | Uniform Metropolis Hastings                 |\n",
    "| $\\mu$           | $(-\\infty,\\infty)$ | Newtonian Monte Carlo (real space proposer) |\n",
    "| $\\sigma$        | $[0,\\infty)$       | Newtonian Monte Carlo (half space proposer) |\n",
    "\n",
    "You can learn more about compositional inference in our framework topics.\n",
    "\n",
    "Normally, this is all you would have to do! However, the HMM model has meaningful\n",
    "structure that we would like to consider when configuring our inference strategy. In\n",
    "particular, the hidden state of each time step, $X_n$, is highly correlated with the\n",
    "state $X_{n-1}$. Thus, we would also like to jointly propose new values for all $X$ —\n",
    "it is very likely that the value of a hidden state $X_n$ becomes invalid and needs to be\n",
    "recomputed after $X_{n-1}$ is updated. Similarly, we would like to update the location\n",
    "$\\mu$ and the scale $\\sigma$ of the hidden states jointly as well. In order to update\n",
    "these variables jointly, we can configure `CompositionalInference` to \"block\" the random\n",
    "variables together.\n",
    "\n",
    "`CompositionalInference` accepts a dictionary that maps families of random variable to\n",
    "the corresponding algorithm, which allow you to override the default inference method\n",
    "for a particular subset of nodes, or group some of them into a block. To define a block,\n",
    "we simply need to pass `CompositionalInference` a *tuple* containing all random variable\n",
    "families that we want to propose jointly as a key. In our case, since we don't want to\n",
    "override the default inference method, we can use `...` instead of providing an\n",
    "inference class for the block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "compositional = bm.CompositionalInference(\n",
    "    {(model.X): ..., (model.Sigma, model.Mu): bm.SingleSiteNewtonianMonteCarlo()}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may notice that we are using what we referred to as \"random variable families\" such\n",
    "as `model.X` as keys, which are essentially functions that generates the random\n",
    "variables, instead of using instantiated random variable like `model.X(0)` and\n",
    "`model.X(1)`. This is because often times the number of random variables is not known\n",
    "ahead of time until an inference starts with some data (you can even have an unbounded\n",
    "number of nodes for some models). By using random variable families in the config, we no\n",
    "longer need to explicitly spell out all instances of the random variables and group them\n",
    "in a huge tuple.\n",
    "\n",
    "The next step is to define the queries and observations. For this particular run, we're\n",
    "interested in inferring $X$, $\\mu$, and $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "queries = (\n",
    "    [model.X(n) for n in range(1, model.N)]\n",
    "    + [model.Mu(k) for k in range(model.K)]\n",
    "    + [model.Sigma(k) for k in range(model.K)]\n",
    ")\n",
    "\n",
    "observations = {\n",
    "    model.X(0): torch.tensor(0.0),\n",
    "    **{model.Y(n): y_obs[n] for n in range(model.N)},\n",
    "    **{model.Theta(k): thetas_obs[k] for k in range(model.K)},\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running inference consists of a few arguments:\n",
    "\n",
    "| Name           | Usage                                                                                                    |\n",
    "| -------------- | -------------------------------------------------------------------------------------------------------- |\n",
    "| `queries`      | List of `@bm.random_variable` targets to fit posterior distributions for.                                |\n",
    "| `observations` | A dictionary of observations, as built above.                                                            |\n",
    "| `num_samples`  | Number of Monte Carlo samples to approximate the posterior distributions for the variables in `queries`. |\n",
    "| `num_chains`   | Number of separate inference runs to use. Multiple chains can help verify that inference ran correctly.  |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6c8e9cf130ff4b3ebe1215367921e4cf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Samples collected:   0%|          | 0/400 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_samples = 400 if not smoke_test else 1\n",
    "num_chains = 1\n",
    "\n",
    "samples = compositional.infer(\n",
    "    queries=queries,\n",
    "    observations=observations,\n",
    "    num_samples=num_samples,\n",
    "    num_chains=num_chains,\n",
    ")\n",
    "\n",
    "samples = samples.get_chain(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X(0) is observed and is not part of query\n",
    "x_samples = torch.stack(\n",
    "    [torch.zeros(num_samples)] + [samples[model.X(n)] for n in range(1, model.N)], dim=1\n",
    ")\n",
    "mu_samples = torch.stack([samples[model.Mu(k)] for k in range(model.K)]).T\n",
    "sigma_samples = torch.stack([samples[model.Sigma(k)] for k in range(model.K)]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "6d7a8437-eddc-434b-a895-325bb761a49f",
    "showInput": false
   },
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "6d7a8437-eddc-434b-a895-325bb761a49f",
    "showInput": false
   },
   "source": [
    "We will look at the values of the samples collected for $X$ and $\\mu$. We will take the\n",
    "mean of samples taken over the last 10% of the chain, and compare these to our synthetic\n",
    "data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEWCAYAAACDoeeyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa40lEQVR4nO3deZRU9Z338fcXGmRfEtRIN5t5ggscdglEJC7RqA+DijwmyigRHeMMqPGEx4HJApnn6NHIRPKIY0aNiiNqPIiacLJgXFBiZhIaoRNEIY+20O2CYjeLgrHl+/xxL6Qaeqnuvrer7s/P6xxOV92q+tWni+pP3/7dW/eauyMiImHqUOgAIiKSHpW8iEjAVPIiIgFTyYuIBEwlLyISMJW8iEjAVPIiIgFTyUu7MzM3s/9R6BwHmNkFZrbNzPaY2ehC5xFJkkpeWszMfm1m/9rA8vPM7G0zKylErjZYBMxx9x7u/lLuDWZWamY1ZjYpZ9mAeNkXmxrUzC4ysxfN7EMzey6d6GBmo81sV+4vTjMba2a1Zja4mcdWmtl2M+ues+zKNPNK+1LJS2ssBf7ezOyQ5ZcCy9y9rgCZ2mIQsLGhG9y9Gvhn4B4z6xIv/g/gPnf/72bGfR9YDNycUM4Gxb+YlgB3W6QTcC/wfXevzGOIjsB1KUaUAlLJS2s8AXwWOOXAAjPrC0wBHjCz8Wb2+3hN8i0zW2JmnRsayMyeM7Mrc65/w8zW5Fw/3syeMrP3zexVM7so57ZzzexlM9ttZtVmNreR5+hgZt81szfitdYHzKy3mR1hZnuISm6Dmf2/Rr7fu4G3gAVmNhM4Dvhucy+Su//W3R8F3mzuvgn4AXAMcBXwL8AeouLPx63AXDPrk040KSSVvLSYu+8FHgUuy1l8EfCKu28APgGuB/oBE4EzgH9q6fPEUwhPAQ8BRwFfB/7dzE6M7/JT4Jvu3hMYDjzTyFDfiP+dBhwL9ACWuPtH7t4jvs9Id/98I9+vA1fG38Ni4B/c/cOWfj9pcvePgCuAW4BvA1e4+/48H74WeA5o8JekZJtKXlprKTA9ZwrjsngZ7l7u7v/l7nXxdMF/AF9uxXNMASrd/b54rJeAx4D/Fd/+MXCimfVy9xp3X9fIODOAH7n7a+6+B5gPfL2F2w7eIFoj3wU834rvpT38GagD/uTur7Twsd8HrjGzI5OPJYWkkpdWcfc1wHvA+Wb2eWA80Ro3ZjbUzFbGG2F3ATcRrdW31CDgi/G0T62Z1RIV9ufi2y8EzgXeMLPVZjaxkXH6E5X0AW8AJcDRLcgyD9gBbKd413j/DVgNlJnZ11vyQHf/M7CS6PuUgGRtLwgpLg8QrcEfB/zG3d+Jl98JvARc7O67zexbwPRGxvgA6JZz/XM5l7cBq939zIYe6O5/BM6LNzTOIZpCGtDAXd8k+oVxwECiNd53GrjvYeLpof8NfBHoDKwxs8fcfUs+j28PZvYVYCpwItEv3PvMbJW7v9+CYRYA64h+WUggtCYvbfEA8BXgH4inamI9iaY19pjZ8cA/NjHGemCamXWLdwG8Iue2lcBQM7vUzDrF/04ysxPMrLOZzTCz3u7+cfx8jc1BPwxcb2ZDzKwH0V8WP8tnLyAz60A09/9Dd3/F3SuA/wvc1cDeRYc+tmM8nVUCdDCzLvEvpETF2y7uAq539/fc/ZdE2zJua8k47v4X4GfAtUlnlMJRyUurxfPtLwLdgZ/n3DQXuATYTbRnys+aGOY24K9Ea9VLgWU54+8GziLa4Pom8DbRhsUj4rtcClTGU0JXE03lNORe4D+J5tJfB/YB1+T3XXId0V8aP8xZ9n+I/uK4ssFH/M2lwF6iv2xOiS/f3dAdzWxg/GGsgfH1GWa2Mef2n5jZTxp5npuINnovy1n2LeAcM2vwr6Am/CvR/6cEwnRmKBGRcGlNXkQkYNrwKtIG8YepGnKOu7/QrmEaEE//vNzIzSe6+9b2zCPtT9M1IiIBK6o1+X79+vngwYMLHUNEJFPKy8vfc/cGP8hWVCU/ePBg1q5dW+gYIiKZYmZvNHabNryKiARMJS8iEjCVvIhIwIpqTr4hH3/8MVVVVezbt6/QUQTo0qULZWVldOqU+KfzRSQFRV/yVVVV9OzZk8GDB9PMoUIkZe7Ojh07qKqqYsiQIYWOIyJ5SLXkzexeomOCb3f34a0ZY9++fSr4ImFmfPazn+Xdd98tdJR2tXl5BdVLVlBSvZW60oGUzpnG0OkjCh2rXWXtNcha3jSlPSd/P3B2WwdRwRePT9v/xeblFWy/YRFWW0PdMWVYbQ3bb1jE5uUVhY7WbrL2GmQtb9pSLXl3f57oZMYimVS9ZAV1vfpCn75Yhw7Qpy91vfpSvWRFoaO1m6y9BlnLm7aC711jZleZ2VozW1us0wA9evRo9j4vvPACw4YNY9SoUezduzfVPAsXLmTRokWpPodESqq34r1611vmvXpTUv3pOeRL1l6DrOVNW8FL3t3vcvdx7j7uyCOze3rJZcuWMX/+fNavX0/Xrl2bvb+7s3///kavS3GoKx2I7dpZb5nt2kld6cACJWp/WXsNspY3bQUv+aRVVMDChTBrVvS1IsFpuOeee45TTz2V6dOnc/zxxzNjxgzcnXvuuYdHH32U733ve8yYEZ234tZbb+Wkk05ixIgRLFiwAIDKykqOO+44LrvsMoYPH84LL7xQ7/q2bdsafBzAjTfeyNChQ5k0aRKvvvpqct+UNKl0zjRKdtVAbQ2+fz/U1lCyq4bSOdMKHa3dZO01yFretBX9LpQtUVEBixZB375QVgY1NdH1uXNhREIb1l966SU2btxI//79Ofnkk/nd737HlVdeyZo1a5gyZQrTp09n1apVbNmyhT/84Q+4O1OnTuX5559n4MCBbNmyhaVLlzJhwgQqKyvrXW/scd27d+eRRx5h/fr11NXVMWbMGMaOHZvMNyRNivbImFtvT42jvnvFp2pPjay9BlnLm7a0d6F8GDgV6GdmVcACd/9pWs+3YkVU8H37RtcPfF2xIrmSHz9+PGVlZQCMGjWKyspKJk2aVO8+q1atYtWqVYwePRqAPXv2sGXLFgYOHMigQYOYMGHCwfvmXm/scbt37+aCCy6gW7fofNdTp05N5puRvAydPuJTWxAHZO01yFreNKVa8u5+cZrjH2rr1mgNPlfv3tHypBxxxBEHL3fs2JG6usPPBe3uzJ8/n29+85v1lldWVtK9e/3TZ+Zeb+xxixcvTiC5iHwaBTUnP3Ag7Ky/vYWdO6Pl7emrX/0q9957L3v2RCcNqq6uZvv27a1+3OTJk3niiSfYu3cvu3fv5he/+EWq+UUkHEHNyU+bFs3BQ7QGv3NnNC9/xRXtm+Oss85i06ZNTJw4EYh2wXzwwQfp2LFjqx43ZswYvva1rzFy5EiOOuooTjrppNS/BxEJQ1Gd/m/cuHF+6ElDNm3axAknnJD3GBUV0Rz81q3RGvy0acnNx0ukpf8nIpIuMyt393EN3RbUmjxEha5SFxGJBDUnLyIi9ankRUQCppIXEQmYSl5EJGAqeRGRgKnk8/DOO+9wySWXcOyxxzJ27FgmTpzI448/3q4ZKisrGT788JNrVVZW8tBDD7VqzMWLF/Phhx8evJ7PIZVFJFtU8s1wd84//3wmT57Ma6+9Rnl5OY888ghVVVWH3behQxykramSby7PoSUvIuEJbj/5pD8N9cwzz9C5c2euvvrqg8sGDRrENddcA8D999/PihUr2LNnD5988gmPP/44s2bN4rXXXqNbt27cddddjBgxgoULF9KjRw/mzp0LwPDhw1m5ciUA55xzDpMmTeLFF1+ktLSUJ598kq5du1JeXs6sWbOA6NOwDZk3bx6bNm1i1KhRzJw5k759+9bL84Mf/IBFixYdfK45c+Ywbtw4du3axZtvvslpp51Gv379ePbZZwH4zne+w8qVK+natStPPvkkRx99dKtfOxEpvLDW5A8ca7impv6xhttwUPmNGzcyZsyYJu+zbt06li9fzurVq1mwYAGjR4+moqKCm266icsuu6zZ59iyZQuzZ89m48aN9OnTh8ceewyAyy+/nNtvv50NGzY0+tibb76ZU045hfXr13P99dcflqcx1157Lf379+fZZ589WPAffPABEyZMYMOGDUyePJm777672ewiUtzCKvncYw136PC3yyuSO7fj7NmzGTlyZL3jx5x55pl85jOfAWDNmjVceumlAJx++uns2LGDXbt2NTnmkCFDGDVqFABjx46lsrKS2tpaamtrmTx5MsDBMfORm6clOnfuzJQpU+rlEJFsC6vkt26NjkyWq43HGh42bBjr1q07eP2OO+7g6aefJvd8tIcePrghJSUl9U7vt2/fvoOX8zl8cUvk5mnqeQ/VqVMnzCyxHCJSeGGVfArHGj799NPZt28fd95558FlTW2sPOWUU1i2bBkQnS6wX79+9OrVi8GDBx/8ZbFu3Tpef/31Jp+3T58+9OnThzVr1gAcHPNQPXv2ZPfu3Y2OM2jQIF5++WU++ugjamtrefrpp/N+rIhkX1glP21aNA9fUwP79//t8rTWn9vRzHjiiSdYvXo1Q4YMYfz48cycOZNbbrmlwfsvXLiQ8vJyRowYwbx581i6dCkAF154Ie+//z7Dhg1jyZIlDB06tNnnvu+++5g9ezajRo2isaOFjhgxgo4dOzJy5Ehuu+22w24fMGAAF110EcOHD+eiiy46eNYpgKuuuoqzzz6b0047LZ+XQkQyKLhDDetYw+nToYZFisun6lDDOtawiMjfhDVdIyIi9WSi5ItpSunTTv8XItlS9CXfpUsXduzYoXIpAu7Ojh076NKlS6GjiEiein5OvqysjKqqqnr7pUvhdOnShbKyskLHEJE8FX3Jd+rUiSFDhhQ6hohIJhX9dI2IiLSeSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQlYST53MrOTgYXAoPgxBri7H5teNBERaau8Sh74KXA9UA58kl4cERFJUr4lv9Pdf5VqEhERSVy+Jf+smd0KrAA+OrDQ3delkkpERBKRb8l/Mf46LmeZA6cnG0dERJKUV8m7+2lN3W5mM919aTKRREQkKUntQnldQuOIiEiCkip5S2gcERFJUFIl7wmNIyIiCdKavIhIwJIq+d8lNI6IiCSoyZI3s9vNrGcDy483s98euO7uc9IIJyIibdPcmvzbwHozuwTAzLqZ2Q+BXwB3pB1ORETapsmSd/cbgTOBGWb2PFAB1AEj3f3xdsgnIiJtkM+c/P74awnQEdjk7h+mF0lERJLS3Jz894DfAg+4+5eAScB5ZrbazE5sj4AiItJ6zR3WoB8w2t13A7h7NTDdzM4BHgNOSDmfiIi0QXNz8tcdKPhDlv8KGJVWKBERSUar95N394+av5eIiBSSzvEqIhKw5ja8Xhd/Pbl94oiISJKaW5O/PP56e9pBREQkec3tXbPJzLYA/c2sIme5Ae7uI9KLJiIibdVkybv7xWb2OeA3wNT2iSQiIklp9vR/7v42MNLMOgND48WvuvvHqSYTEZE2y+scr2b2ZeABoJJoqmZAfF7X51PMJiIibZRXyQM/As5y91cBzGwo8DAwNq1gIiLSdvnuJ9/pQMEDuPtmoFM6kUREJCn5rsmvNbN7gAfj6zOAtelEEhGRpORb8v8IzAauja+/APx7KolERCQxeZV8fJyaH8X/REQkI3TsGhGRgKnkRUQC1qKSN7NuaQUREZHk5VXyZvYlM3sZeCW+PtLMtOFVRKTI5bsmfxvwVWAHgLtvACanFUpERJKR93SNu287ZNEnCWcREZGE5buf/DYz+xLgZtYJuA7YlF4sERFJQr5r8lcTfRiqFKgmOon37JQyiYhIQvL9MNR7RIcyEBGRDMn3UMP3AX7ocneflXgiERFJTL5z8itzLncBLgDeTD6OiIgkKd/pmsdyr5vZw8CaVBKJiEhiWntYgy8ARyUZREREkpfvnPxuojl5i7++DfxzirlERCQB+U7X9Ew7iIiIJK/JkjezMU3d7u7rko0jIiJJam5N/t+auM2B0xPMIiIiCWuy5N39tPYKIiIiyct3P3nMbDhwItF+8gC4+wNphBIRkWTku3fNAuBUopL/JXAO0X7yKnkRkSKW737y04EzgLfd/XJgJNA7tVQiIpKIfEt+r7vvB+rMrBewHRiQXiwREUlCvnPya82sD3A3UA7sAX6fVigREUlGc/vJ3wE85O7/FC/6iZn9Gujl7hWppxMRkTZpbk1+M7DIzI4BHgUedveX0o8lIiJJaHJO3t1/7O4TgS8TncT7XjN7xcwWmNnQdkkoIiKtlteGV3d/w91vcffRwMXA+egcryIiRS+vkjezEjP7OzNbBvwKeBWYlmoyERFps+Y2vJ5JtOZ+LvAH4BHgKnf/oB2yiYhIGzW34XU+8BDwbXevaYc8IiKSoOYOUKajTIqIZFhrT/8nIiIZoJIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQmYSl5EJGAqeRGRgKnkRUQCppIXEQlYSdpPYGZnAz8GOgL3uPvNST/H5uUVVC9ZQUn1VupKB1I6ZxpDp48oujE1bvaypilLr21aspQ3q+/bVNfkzawjcAdwDnAicLGZnZjkc2xeXsH2GxZhtTXUHVOG1daw/YZFbF5eUVRjatzsZU1Tll7btGQpb5bft2lP14wH/uLur7n7X4FHgPOSfILqJSuo69UX+vTFOnSAPn2p69WX6iUrimpMjZu9rGnK0mublizlzfL7Nu2SLwW25VyvipcdZGZXmdlaM1v77rvvtvgJSqq34r1611vmvXpTUr21FXHTG1PjpjdmmuOmJUuvbVqylDfL79uCb3h197vcfZy7jzvyyCNb/Pi60oHYrp31ltmundSVDmx1pjTG1LjpjZnmuGnJ0mublizlzfL7Nu2SrwYG5Fwvi5clpnTONEp21UBtDb5/P9TWULKrhtI504pqTI2bvaxpytJrm5Ys5c3y+9bcPbHBDhvcrATYDJxBVO5/BC5x940N3X/cuHG+du3aFj9PlvZS0LjZypqmLL22aclS3mJ+35pZubuPa/C2NEs+fvJzgcVEu1De6+43Nnbf1pa8iMinWVMln/p+8u7+S+CXaT+PiIgcruAbXkVEJD0qeRGRgKnkRUQCppIXEQlY6nvXtISZvQu80YYh+gHvJRQnbVnKCtnKm6WskK28WcoK2crblqyD3L3BT5MWVcm3lZmtbWw3omKTpayQrbxZygrZypulrJCtvGll1XSNiEjAVPIiIgELreTvKnSAFshSVshW3ixlhWzlzVJWyFbeVLIGNScvIiL1hbYmLyIiOVTyIiIBC6LkzexsM3vVzP5iZvMKnacpZjbAzJ41s5fNbKOZXVfoTM0xs45m9pKZrSx0luaYWR8zW25mr5jZJjObWOhMjTGz6+P3wJ/N7GEz61LoTLnM7F4z225mf85Z9hkze8rMtsRf+xYy4wGNZL01fh9UmNnjZtangBHraShvzm3fNjM3s35JPFfmS749ThaesDrg2+5+IjABmF3keQGuAzYVOkSefgz82t2PB0ZSpLnNrBS4Fhjn7sOJDsX99cKmOsz9wNmHLJsHPO3uXwCejq8Xg/s5POtTwHB3H0F0Xov57R2qCfdzeF7MbABwFpDY+f8yX/K0w8nCk+Tub7n7uvjybqISKm36UYVjZmXA/wTuKXSW5phZb2Ay8FMAd/+ru9cWNFTTSoCu8cl1ugFvFjhPPe7+PPD+IYvPA5bGl5cC57dnpsY0lNXdV7l7XXz1v4jOTFcUGnltAW4DbgAS2yMmhJJv9mThxcrMBgOjgf8ucJSmLCZ60+0vcI58DAHeBe6Lp5fuMbPuhQ7VEHevBhYRrbG9Bex091WFTZWXo939rfjy28DRhQzTArOAXxU6RFPM7Dyg2t03JDluCCWfSWbWA3gM+Ja77yp0noaY2RRgu7uXFzpLnkqAMcCd7j4a+IDimU6oJ57LPo/oF1N/oLuZ/X1hU7WMR/tfF/0+2Gb2HaJp0mWFztIYM+sG/Avw/aTHDqHkUz9ZeNLMrBNRwS9z9xWFztOEk4GpZlZJNA12upk9WNhITaoCqtz9wF9Gy4lKvxh9BXjd3d9194+BFcCXCpwpH++Y2TEA8dftBc7TJDP7BjAFmOHF/aGgzxP9wt8Q/7yVAevM7HNtHTiEkv8j8AUzG2JmnYk2Xv28wJkaZWZGNGe8yd1/VOg8TXH3+e5e5u6DiV7XZ9y9aNc23f1tYJuZHRcvOgN4uYCRmrIVmGBm3eL3xBkU6UbiQ/wcmBlfngk8WcAsTTKzs4mmGqe6+4eFztMUd/+Tux/l7oPjn7cqYEz8nm6TzJd8vGFlDvAboh+SR919Y2FTNelk4FKiteL18b9zCx0qINcAy8ysAhgF3FTYOA2L/9pYDqwD/kT0s1hUH8E3s4eB3wPHmVmVmV0B3AycaWZbiP4aubmQGQ9oJOsSoCfwVPxz9pOChszRSN50nqu4/4IREZG2yPyavIiINE4lLyISMJW8iEjAVPIiIgFTyYuIBEwlLyISMJW8iEjAVPIizTCzwfGx6e+Oj/++ysy6FjqXSD5U8iL5+QJwh7sPA2qBCwsbRyQ/KnmR/Lzu7uvjy+XA4MJFEcmfSl4kPx/lXP6E6LDGIkVPJS8iEjCVvIhIwHQUShGRgGlNXkQkYCp5EZGAqeRFRAKmkhcRCZhKXkQkYCp5EZGAqeRFRAL2/wEF7v2CfB5nBwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEWCAYAAACAOivfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZdUlEQVR4nO3dfbiVdZ3v8fdXBEFBoEQbRR7yhCVcG5QtgyNSOqNZhyzJzDK1tEMPWh2v4cxonY40Z/LUxJRdZc4hw4fxaTyJWp5KHJ/pYQwQKCTDY1vcPiQhIBiYyPf8sRa0QdisDWutm73v9+u69rXXfa+1fr/vgutan33/fvf9uyMzkSSV0z5FFyBJKo4hIEklZghIUokZApJUYoaAJJWYISBJJWYISFKJGQJqmojIiPhPRdexRUScHhFPR8T6iDi66HqkIhgCqllE/CQi/mEH+98bEc9HxL5F1LUHZgIXZWb/zHy0iAIi4sSIuD8i1kZEWxE11FNEXBsR/9hhe3REPBcR04usSztnCKgrrgM+EhGx3f5zgBszc1MBNe2J4cDSgmt4GZgN/LeC66i76tHV/cA/ZubMouvRjhkC6oo7gDcCJ2zZERGDgSnA9RExISJ+HhFrqn/9fTsi+uyooYh4ICI+3mH7oxExr8P2WyPinoh4MSIej4gzOzz37oh4LCLWRcQzO/srMyL2iYj/HhFPRcQLEXF9RAyMiP0iYj3QC1gcEf9vJ+/PiPh0RCyv9vU/I+KIiPhZRLwUEbdu+Xzb19/h/Z0Of2XmI5n5r8CTnb2uHprxeTq8dgJwD/D5zLyy/p9G9WIIqGaZuQG4FTi3w+4zgd9k5mLgNeBi4CDgOOCvgU93tZ+IOIDKF8hNwMHAWcB3IuKo6ku+B3wiMwcAY4D7dtLUR6s/JwJvBvoD387MVzKzf/U1YzPziE7KeScwHpgI/B0wC/gIcHi17w919fMVrBmfZwLwE+DizLy6Du2pgQwBddV1wBkR0be6fW51H5m5IDN/kZmbMrMN+N/A23ejjylAW2ZeU23rUeA24APV518FjoqIAzNzdWYu3Ek7ZwNfz8wnM3M9cClwVhfnLv4pM1/KzKXAr4G51fbWAj8GutuEcjM+z0RgS3vayxkC6pLMnAf8AXhfRBxB5a++mwAiYlRE3FWdJH4JuJzKUUFXDQf+sjqstCYi1lD5Qn9T9fn3A+8GnoqIByPiuJ20cyjwVIftp4B9gUO6UMvvOzzesIPt/nQvzfg8VwLzgXuqw4XaixkC2h3XUzkC+Ahwd2Zu+SK5CvgN8JbMPBD4PLD9JPIWLwP7d9h+U4fHTwMPZuagDj/9M/NTAJn5y8x8L5WhojuoDFHtyLNUAmWLYcAmtv3iq5dtPk9EvKmT13YHe/J5XgM+DKwA7o6IA+tcm+rIENDuuB74G+C/UB0KqhoAvASsj4i3Ap/qpI1FwNSI2L862XhBh+fuAkZFxDkR0bv6c2xEvC0i+kTE2RExMDNfrfa3eSd93AxcHBEjI6I/lSOTf2vQWUyLgdERMa46VDajljdVJ6/7Ar0rm9F3Z5PpTbZbn2eL6v/NB6gcNf6oOs+jvZAhoC6rjvf/DDgA+EGHp6ZT+QtwHfBd4N86aeYbwJ+o/FV+HXBjh/bXAadQmRB+Fnge+CqwX/Ul5wBt1SGnT1IZKtqR2cC/Ag8BvwM2Ap+p7VN2TWb+FvgH4N+B5cC8zt+x1WQqwzA/onKksgGYu7MXR+XCthOqj0+onuW05bnPR0RdxuH34PN0bONPwFQq/+4/jIh+9ahN9RXeWUySyssjAUkqMUNAarCIWFodxtn+Z2fDWHu1nvZ5ys7hIEkqse624BcHHXRQjhgxougyJKlbWbBgwR8yc8j2+7tdCIwYMYL58+cXXYYkdSsR8dSO9jsnIEklZghIUokZApJUYt1uTmBHXn31Vdrb29m4cWPRpZRe3759GTp0KL179y66FEk16BEh0N7ezoABAxgxYgTxupteqVkyk1WrVtHe3s7IkSOLLkdSDXpECGzcuNEA2AtEBG984xtZuXJl0aVIPceSJTBnDqxYAcOGwdSp0NJSt+Z7zJyAAbB38P9BqqMlS2DmTFi9GoYOrfyeObOyv056xJGAJPVIc+aw8rXBLF08mLVrYeDAwYw+DIbMmVO3o4EecyRQtP79d31DpocffpjRo0czbtw4NmzY0NB6ZsyYwcyZMxvah6TGenHRCub9aiAbNsCBB8KGDTDvVwN5cdGKuvVhCDTRjTfeyKWXXsqiRYvo12/XS6tnJps3b97ptqSebfGaYQzeZy39+kEE9OsHg/dZy+I1w+rWRylDYMkSmDEDzj+/8ruOw2s88MADvOMd7+CMM87grW99K2effTaZydVXX82tt97KF7/4Rc4+u7LY4te+9jWOPfZYWlpauOyyywBoa2vjyCOP5Nxzz2XMmDE8/PDD22w//fTTO3wfwJe//GVGjRrFpEmTePzxx+v3oSQV4r5BUxm4eTV9N6yG3EzfDasZuHk19w2aWrc+SjcnsGWeZfDgbedZpk+v34T7o48+ytKlSzn00EM5/vjj+elPf8rHP/5x5s2bx5QpUzjjjDOYO3cuy5cv55FHHiEzOe2003jooYcYNmwYy5cv57rrrmPixIm0tbVts72z9x1wwAHccsstLFq0iE2bNnHMMccwfvz4+nwgSYXoNa6Fu/efzl8+O4dBa1ewZuAwHjjiAnqNqt/ZQaULgTlzKgEweHBle8vvOs6zMGHCBIYOHQrAuHHjaGtrY9KkSdu8Zu7cucydO5ejjz4agPXr17N8+XKGDRvG8OHDmThx4tbXdtze2fvWrVvH6aefzv77V+4Nftppp9Xnw0gqzNSpMHNmC8+NbWHgQFi7tvKH6/T6HQiULwRWrKgcAXQ0cGBlf73st99+Wx/36tWLTZtef1/zzOTSSy/lE5/4xDb729raOOCAbe/J3XF7Z++74oor6lC5pL1JS0tllKLjZQIXXFDXywTKNycwbFglTTtau7ayv5ne+c53Mnv2bNavr9wn/JlnnuGFF17Y7fdNnjyZO+64gw0bNrBu3Tp++MMfNrR+Sc3R0lKZu5w9u/K7ngEAJTwSqBxeVR53PLy64ILm1nHKKaewbNkyjjvuOKByiukNN9xAr169dut9xxxzDB/84AcZO3YsBx98MMcee2zDP4Ok7q/b3V6ytbU1t7+pzLJly3jb295WcxsNvgq79Lr6/yGp8SJiQWa2br+/dEcCUPnC90tfkko4JyBJ+jNDQJJKzBCQpBIzBCSpxPaKEIiIXhHxaETcVXQtklQme0UIAJ8DlhVdxJ74/e9/z4c//GHe/OY3M378eI477jhuv/32ptbQ1tbGmDFjdrj/pptu2q02r7jiCv74xz9u3a5lyWxJ3UfhIRARQ4H/DFxddC27KzN53/vex+TJk3nyySdZsGABt9xyC+3t7a977Y6WkGi0zkJgV/VsHwKSepa94TqBK4C/Awbs7AURMQ2YBjCsHus71Plqsfvuu48+ffrwyU9+cuu+4cOH85nPfAaAa6+9ljlz5rB+/Xpee+01br/9ds4//3yefPJJ9t9/f2bNmkVLSwszZsygf//+TJ8+HYAxY8Zw112VEbJ3vetdTJo0iZ/97Gccdthh3HnnnfTr148FCxZw/vnnA5WriXfkkksuYdmyZYwbN47zzjuPwYMHb1PPl770JWbOnLm1r4suuojW1lZeeuklnn32WU488UQOOugg7r//fgC+8IUvcNddd9GvXz/uvPNODjnkkN3+t5NUrEKPBCJiCvBCZi7o7HWZOSszWzOzdciQIXvWaQPu2bl06VKOOeaYTl+zcOFCvv/97/Pggw9y2WWXcfTRR7NkyRIuv/xyzj333F32sXz5ci688EKWLl3KoEGDuO222wD42Mc+xre+9S0WL1680/d+5Stf4YQTTmDRokVcfPHFr6tnZz772c9y6KGHcv/9928NgJdffpmJEyeyePFiJk+ezHe/+91d1i5p71X0cNDxwGkR0QbcApwUETc0tMeOa0nvs8+fH8+ZU7cuLrzwQsaOHbvN+j0nn3wyb3jDGwCYN28e55xzDgAnnXQSq1at4qWXXuq0zZEjRzJu3DgAxo8fT1tbG2vWrGHNmjVMnjwZYGubtehYT1f06dOHKVOmbFOHpO6r0BDIzEszc2hmjgDOAu7LzI80tNMVKyorx3W0h2tJjx49moULF27dvvLKK7n33ntZuXLl1n3bLw+9I/vuu+82t4/cuHHj1se1LE/dFR3r6azf7fXu3ZuIqFsdkopV9JFA8zVgLemTTjqJjRs3ctVVV23d19lk6gknnMCNN94IVG5HedBBB3HggQcyYsSIrWGycOFCfve733Xa76BBgxg0aBDz5s0D2Nrm9gYMGMC6det22s7w4cN57LHHeOWVV1izZg333ntvze+V1L3tNSGQmQ9k5pSGdzR1amUeYPVq2Lz5z4+n7v6teiKCO+64gwcffJCRI0cyYcIEzjvvPL761a/u8PUzZsxgwYIFtLS0cMkll3DdddcB8P73v58XX3yR0aNH8+1vf5tRo0btsu9rrrmGCy+8kHHjxrGzFWFbWlro1asXY8eO5Rvf+Mbrnj/88MM588wzGTNmDGeeeebWu5YBTJs2jVNPPZUTTzyxln8KSd1MKZeSdi3pxnIpaWnv41LSHbmWtCQBe9FwkCSp+XpMCHS3Ya2eyv8HqXvpESHQt29fVq1a5RdQwTKTVatW0bdv36JLkVSjHjEnMHToUNrb27c5L1/F6Nu3L0OHDi26DEk16hEh0Lt3b0aOHFl0GZLU7fSI4SBJ0u4xBCSpxAwBSSoxQ0CSSswQkKQSMwQkqcQMAUkqMUNAkkrMEJCkEjMEJKnEDAFJKjFDQJJKzBCQpBIzBCSpxAwBSSoxQ0CSSswQkKQSMwQkqcQMAUkqMUNAkkrMEJCkEis0BCKib0Q8EhGLI2JpRHypyHokqWz2Lbj/V4CTMnN9RPQG5kXEjzPzFwXXJUmlUGgIZGYC66ubvas/WVxFklQuhc8JRESviFgEvADck5n/sYPXTIuI+RExf+XKlU2vUZJ6qsJDIDNfy8xxwFBgQkSM2cFrZmVma2a2DhkypOk1SlJPVXgIbJGZa4D7gVMLLkWSSqPos4OGRMSg6uN+wMnAb4qsSZLKpOizg/4CuC4ielEJpFsz866Ca5Kk0ij67KAlwNFF1iBJZbbXzAlIkprPEJCkEjMEJKnEDAFJKjFDQJJKrKYQiIjYwb796l+OJKmZaj0S+F7HjYjoD/yo/uVIkpqp1hBoj4jvAETEYGAucEPDqpIkNUVNIZCZ/wNYHxH/QiUA/jkzr2loZZKkhuv0iuGImNph8z+ALwKPABkRUzNzTiOLkyQ11q6WjXjPdtuPUrnxy3uo3PzFEJCkbqzTEMjMj9XSSERcmpn/qz4lSZKapV7XCXygTu1IkpqoXiHwuusIJEl7v3qFgDeHl6RuyCMBSSqxeoXA/6lTO5KkJqrpzmIRcQ07GPLJzPOrvy+vc12SpCao9faSHe/72xc4HXi2/uVIkpqpphDIzNs6bkfEzcC8hlQkSWqa3Z0TeAtwcD0LkSQ1X61zAuuozAlE9ffzwN83sC5JUhPUOhw0oNGFSJKar9aJYSKiBRjR8T2uIipJ3Vutw0GzgRZgKbC5uttVRCWpm6v1SGBiZh7V0EokSU1X69lBP48IQ0CSephajwSupxIEzwOvUD1LKDNb9qTziDi82vYhVIaXZmXmN/ekTUlS7WoNge8B5wC/4s9zAvWwCfjbzFwYEQOABRFxT2Y+Vsc+JEk7UWsIrMzMH9S788x8Dniu+nhdRCwDDgMMAUlqglpD4NGIuAn4IZXhIKC+p4hGxAjgaCo3tN/+uWnANIBhw4bVq0tJKr1aQ6AflS//Uzrsq9spohHRH7gN+K+Z+dL2z2fmLGAWQGtrqzewkaQ6qfWK4ZpuOL87IqI3lQC40YvPJKm56nVTmd0SEUFl0nlZZn69yFokqYwKDQHgeCpnHZ0UEYuqP+8uuCZJKo1Oh4Mi4nOZ+c2IOD4zf1rvzjNzHt6fWJIKs6sjgS1zAd9qdCGSpObb1cTwsohYDhwaEUs67K/LFcOSpGJ1GgKZ+aGIeBNwN3Bac0qSJDXLLk8RzczngbER0QcYVd39eGa+2tDKJEkNV+v9BN5OZaG3NipDQYdHxHmZ+VADa5MkNVitVwx/HTglMx8HiIhRwM3A+EYVJklqvFqvE+i9JQAAMvO3QO/GlCRJapZajwTmR8TVwA3V7bOB+Y0pSZLULLWGwKeAC4HPVrcfBr7TkIokSU1T6wJyr1CZF3B9H0nqQYpeO0iSVCBDQJJKrEshEBH7N6oQSVLz1RQCEfFXEfEY8Jvq9tiIcGJYkrq5Wo8EvgG8E1gFkJmLgcmNKkqS1Bw1Dwdl5tPb7XqtzrVIkpqs1usEno6IvwKyek/gzwHLGleWJKkZaj0S+CSVi8UOA54BxlW3JUndWK0Xi/2BylIRkqQepNalpK8Bcvv9mXl+3SuSJDVNrXMCd3V43Bc4HXi2/uVIkpqp1uGg2zpuR8TNwLyGVCRJaprdXTbiLcDB9SxEktR8tc4JrKMyJxDV388Df9/AuiRJTVDrcNCARhciSWq+TkMgIo7p7PnMXFjfciRJzbSrI4F/7uS5BE6qYy2SpCbrNAQy88RmFSJJar5arxMgIsYAR1G5TgCAzLx+TwuIiNnAFOCFzByzp+1JkmpX6/0ELgO+Vf05Efgn4LQ61XAtcGqd2pIkdUGt1wmcAfw18HxmfgwYCwysRwGZ+RDwYj3akiR1Ta0hsCEzNwObIuJA4AXg8MaVta2ImBYR8yNi/sqVK5vVrST1eLWGwPyIGAR8F1gALAR+3qiitpeZszKzNTNbhwwZ0qxuJanH29V1AlcCN2Xmp6u7/iUifgIcmJlLGl6dJKmhdnV20G+BmRHxF8CtwM2Z+Wjjy5IkNUOnw0GZ+c3MPA54O5WbzM+OiN9ExGURMaoeBVRXJP05cGREtEfEBfVoV5K0a5H5unvFdP6GiKOB2UBLZvZqSFWdaG1tzfnz5ze7W0nq1iJiQWa2br+/1usE9o2I90TEjcCPgceBqXWuUZLUZLuaGD4Z+BDwbuAR4BZgWma+3ITaJEkNtquJ4UuBm4C/zczVTahHktREu1pAzlVCJakH293bS0qSegBDQJJKzBCQpBIzBCSpxAwBSSoxQ0CSSswQkKQSMwQkqcQMAUkqMUNAkkrMEJCkEjMEJKnEDAFJKjFDQJJKzBCQpBIzBCSpxAwBSSoxQ0CSSswQkKQSMwQkqcQMAUkqMUNAkkrMEJCkEis8BCLi1Ih4PCKeiIhLiq5Hksqk0BCIiF7AlcC7gKOAD0XEUUXWJEllUvSRwATgicx8MjP/BNwCvLfgmiSpNIoOgcOApztst1f3bSMipkXE/IiYv3LlyqYVJ0k9XdEhUJPMnJWZrZnZOmTIkKLLkaQeo+gQeAY4vMP20Oo+SVITFB0CvwTeEhEjI6IPcBbwg4JrkqTS2LfIzjNzU0RcBNwN9AJmZ+bSImuSpDIpNAQAMvNHwI+KrkOSyqjo4SBJUoEMAUkqMUNAkkrMEJCkEjMEJKnEDAFJKjFDQJJKzBCQpBIzBCSpxAwBSSoxQ0CSSswQkKQSMwQkqcQMAUkqMUNAkkrMEJCkEjMEJKnEDAFJKjFDQJJKzBCQpBIzBCSpxAwBSSoxQ0CSSswQkKQSMwQkqcT2LbqApliyBObMgRUrYNgwmDoVWlqKrkqSCtfzQ2DJElZ/YSaPvzCY378ylEOWrubIBTMZ/OXpBoGk0itsOCgiPhARSyNic0S0Nqqf56+awyNPDGY1gzlw4D6sZjCPPDGY56+a06guJanbKHJO4NfAVOChRnby3C9WsHnAQPr1gwjo1w82DxjIc79Y0chuJalbKCwEMnNZZj7e6H6eymEMYu02+waxlqdyWKO7lqS9Xo8/O+jZiVPptW41fTeshtxM3w2r6bVuNc9OnFp0aZJUuIaGQET8e0T8egc/7+1iO9MiYn5EzF+5cmWXapj06RbuOGJ6ZU5gbTurGcwdR0xn0qedFJakhp4dlJl/U6d2ZgGzAFpbW7Mr721pAS5vYc6clq1niJ7lGaKSBJThFFEqX/h+6UvS6xV5iujpEdEOHAf834i4u6haJKmsCjsSyMzbgduL6l+SVIKzgyRJO2cISFKJGQKSVGKR2aUzLgsXESuBp3bz7QcBf6hjOZLULHv6/TU8M4dsv7PbhcCeiIj5mdmwxeokqVEa9f3lcJAklZghIEklVrYQmFV0AZK0mxry/VWqOQFJ0rbKdiQgSerAEJCkEitNCETEqRHxeEQ8ERGXFF2PJNUiImZHxAsR8etGtF+KEIiIXsCVwLuAo4APRcRRxVYlSTW5Fji1UY2XIgSACcATmflkZv4JuAXo0t3NJKkImfkQ8GKj2i9LCBwGPN1hu726T5JKrSwhIEnagbKEwDPA4R22h1b3SVKplSUEfgm8JSJGRkQf4CzgBwXXJEmFK0UIZOYm4CLgbmAZcGtmLi22KknatYi4Gfg5cGREtEfEBXVt32UjJKm8SnEkIEnaMUNAkkrMEJCkEjMEJKnEDAFJKjFDQNpDETGiUSs8So1mCEhSiRkCUh1FxJsj4tGIOLboWqRa7Ft0AVJPERFHUlmm/KOZubjoeqRaGAJSfQwB7gSmZuZjRRcj1crhIKk+1gIrgElFFyJ1hUcCUn38CTgduDsi1mfmTUUXJNXCEJDqJDNfjogpwD3VIHC5cu31XEVUkkrMOQFJKjFDQJJKzBCQpBIzBCSpxAwBSSoxQ0CSSswQkKQS+/92C6FBj6MoGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tail_len = num_samples // 10\n",
    "\n",
    "xs_tail_plot = x_samples[-tail_len:].mean(0)\n",
    "plt.scatter(\n",
    "    range(len(xs_tail_plot)),\n",
    "    xs_tail_plot,\n",
    "    alpha=0.5,\n",
    "    label=\"Inferred\",\n",
    "    c=\"b\",\n",
    ")\n",
    "plt.scatter(range(len(xs_obs)), xs_obs, alpha=0.5, label=\"Ground truth\", c=\"r\")\n",
    "plt.yticks(range(K))\n",
    "plt.title(\"Values of X_1 ... X_N\")\n",
    "plt.xlabel(\"n\")\n",
    "plt.ylabel(\"Value of X_n\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "mus_tail_plot = mu_samples[-tail_len:].mean(0)\n",
    "plt.scatter(range(K), mus_tail_plot, alpha=0.5, label=\"Inferred\", c=\"b\")\n",
    "plt.scatter(range(K), mus_obs, alpha=0.5, label=\"Ground truth\", c=\"r\")\n",
    "plt.xticks(range(K))\n",
    "plt.title(\"Values of mu_1 ... mu_K\")\n",
    "plt.xlabel(\"k\")\n",
    "plt.ylabel(\"Value of mu_k\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "e9beca04-43f5-4271-bf4d-90859b2a10d3",
    "showInput": false
   },
   "source": [
    "These plots indicate that inference seems to be recovering hidden states well, and is\n",
    "computing reasonably accurate values for $\\mu$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "4ea45674-027b-4307-b1af-5d1a10cad813",
    "showInput": false
   },
   "source": [
    "## Posterior likelihood checks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "4ea45674-027b-4307-b1af-5d1a10cad813",
    "showInput": false
   },
   "source": [
    "One way to evaluate posterior samples is computing likelihoods of our posterior samples,\n",
    "and comparing these to the likelihood of the underlying synthetic data. Formally, we can\n",
    "compute the joint likelihood of posterior samples with the observations used to generate\n",
    "the samples. And similarly, we can compute the joint likelihood of the observations with\n",
    "the underlying synthetic data which was generated at the same time as the observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(xs, ys, thetas, mus, sigmas, N):\n",
    "    \"\"\"Returns the log likelihood of the HMM model conditioned on the data\"\"\"\n",
    "    result = 0\n",
    "    # transition probabilities\n",
    "    for n in range(1, N):\n",
    "        result += torch.log(thetas[xs[n - 1], xs[n]])\n",
    "    # emission probabilities\n",
    "    for n in range(N):\n",
    "        result += dist.Normal(mus[xs[n]], sigmas[xs[n]]).log_prob(ys[n])\n",
    "    return result\n",
    "\n",
    "\n",
    "# computes the log likelihood of the HMM model per iteration\n",
    "ppcs = [\n",
    "    log_likelihood(x, y_obs, thetas_obs, mu, sigma, N)\n",
    "    for x, mu, sigma in zip(x_samples.int(), mu_samples, sigma_samples)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAFlCAYAAADYhD9JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+mElEQVR4nO3deXhV1b3/8c83MxDmGcIQ5hmEgCCoyKA4AXq1zji1VKu0UjtovVXrhTrWqVVbVKq2DtefAyJaB0RUFEEmkUEQkFHGhDEk5Jyc9fsjJ7kBkpBE9tkh+/16Hh6y197n7G/WswOfrLP22uacEwAAAABvxPldAAAAAFCdEbgBAAAADxG4AQAAAA8RuAEAAAAPEbgBAAAADxG4AQAAAA8l+F2A1xo1auTatm3rdxkAAACoxhYuXLjLOde4pH3VPnC3bdtWCxYs8LsMAAAAVGNmtqG0fUwpAQAAADxE4AYAAAA8ROAGAAAAPETgBgAAADxE4AYAAAA8ROAGAAAAPHTCBW4zG2Vmq8xsjZnd5nc9AAAAQFlOqMBtZvGSnpB0tqRuki4zs27+VgUAAACU7oQK3JIGSFrjnFvnnMuT9IqkMT7XBAAAAJTqRAvcLSVtKra9Odp2GDMbb2YLzGzBzp07Y1YcAAAAcKQTLXCXi3NuinMuwzmX0bhxiY+0BwAAAGIiwe8CKmiLpFbFttOibQCOA+ecdh7cqSa1mvhdCk4gefl5+nTDpwrlhyRJDWo00MlpJ/tclfciLqI4q5bjVgCOsxMtcH8lqaOZpasgaF8q6XJ/SyrBLbdIS5b4XQUCwEnKzjugWkm1ZLLD2q3UV5Us4pxWZ67WtgPblNC0lxrUqF/uGg6GshWO5KtOcu3D6ii0O3e3DoZylBiXqMa1Gle4NlRtO/ZtUULWmqL/UHIkHWjRT6lJqUXH5IYPafnOZeraqJtqJtbwpc4fa9+h/dq8b5M6N+qizIOZ+i5rtXo06am6yXWO+VonJ+cKvjKLO+bPQMRFtCN7p/Ly8w5rT4xPVLPUpiX+nEUKTqA4K/3dtx3Yri37Nx/WZjK1qddGDWs0LMf3cex/W5ykTXs3KSsnSzIpKS5JifH/FzdSk2qreWqzY54LqLQ+faRHH/W7isOcUIHbORc2s5slvS8pXtJU59xyn8sCYspJ2p2zW/sO7dOOgzt0MO+gOjfqpOapzSVJG/du0g/7t+ik5icpOT656HU7D+7Smqw1ci5S4vtGXEThSL7i4+K1bvda1a+RcdR/rHn5ecp3+aqRUEMb927Uhr0b5ZxTJPqeSQlJal2nlVrUbln0n35u+JCWbl8aDRtSw+wGalO3rQozQc3EWopnlPCElpWbpRqJKeraqKvyXURfb/9aOw/uOixw7z20R/sPHdCmfRvVuWHnw16fG85VckKKp7+IHQwd1OrM79ShQQelJtWq1Htsz96uHdk7lZqUqq0HtimUH9ayHcvUr3lfpSSklPq6lbu+1fYD24u2U5NqqXezPkqMO/y/4AN52dq4d4PyIiEdzMtWXvQTgyNtO7BN3Rt3U1J8UlFbxEW0aOtiJcYnqFfT3iX2Zb7L19rdaxQfl6CaiTWL2nNCOVq+Y7l6Ne2teil1S/0+dufu0bIdy5Rer61a1kmTqeDfo4jLV34kX7tz9+hQOFf7Du3TroOZSk2qpbi4eO3P26dQJFxQZyRfsm0EbgTOCRW4Jck5966kd/2uo0xV7LcqVC///vpfGjdtnCSpb/O+2rR3k05r01Ov/eQ1SdIFU/pp0dZ1ymgR1uyrZ6pWUi1t2bdF3Z/srrQ6HTS41eBS33tUh1HKCefoijeu0B+GDNZJzU/SiHYjVC+lnlZnrtbpz52ucCSs1TevVr+/dVJanZ4a1naYejTpoZSEFD27+Fl99P1HOqlZHX1w1QdqVLORfvefX+qpBQu0aPwifbLhE018f6LCkUVF57yq11V64YIXvO00eCYvP0+n3N9A43pfqyfPfVKSdNdzp2tP7h59fcPsouOe/Pgu3fPpPUqO362NE18tmra0Yc8GpT+Wrl+d/Cs9MuoRT2oMR8I6Y+pgzd+yR6e0StKcaz+WlTEKXJornh6gr37YonjbqHyXr8nDJuv+z+9X/xapmjluZomvWbp9qXr/vbfGdhmr/i3661D4kO77/D71ax5Ro5p1tHDrwqJjt+7fqropddWjSQ81rtlYEwZM0ICWAw57vze/fVPnzvi5Gtf8Qe9f+b46Nyr45eUPM2/T/Z9/Jkl6/8o/6KstX2naqmlKjk/W7UNu17mdztXDnz+o382co7nXf6ZeaQOL3jPzYKYGTx2szJzN2vLreYcF+eLueW+iHp33taS1+sOQSzR5+GRd9Op/6Y2Vbxx2XGJcoiYNu1+/PeW3R/XzA58/oN/P/L323z7jsF/IgOruhAvcQJDtzN6pie9P1KC0Qfrgqg+UmpSq69+6Xm9++6byI/nKzMnUoq2LNDx9uGZ9P0tpj6RpZLuRWpW5SqFISNMunaYODTqUeY6Ii+hv8/+mP8/5syQpJSFFfZv31erM1crLz9O+Q/t00f+7SLsO7tK/L/i3zupwVtFrL+1xqV5f+bqufONKnfPiObp/xP16ZtEzuqrXVerZtKd6Nu2pEe1G6LvM7yRJkz6bpGU7lnnXYfDc/C3zlR3K1oh2I4raxnYeq19/8Gut271O7eq3kySt2b1GtZNqa3/efj311VO6a+hdkqTVmavl5PTovEdlZkXHH0vf5n01KG1QicE5HAnriflPaNG2gl/stu7fqvlb5uuibhfptRWvafJnkzUobZCGth2q+Lj4cp0vJ5SjxdsWa1SHUXpvzXvq1ribbhtym2om1tTE9yfqo3UfaXi74Ue97p5P7lGd5DqaOnqq6kenaXVp1EWXv3G5mqU201ntz1JiXKIkKa1Omn558i+LjivJlb2uVLfG3XT2i2erzz/6qEGNBkXf47je4zTr+1m6/PXLlZmTqYFpA7Uje4fOe/k8jWg3Qgt+WKAz25+pgcXCtiQ1rNlQ9w6/Vxe+eqHmbZ6nU9ucWuK5526eq8GtBis5IVlvrXpL/zPsf/TB2g80tO1Qjek8RoNbDVb3Jt2VEJdQamhvFh3Z3n5gu1IbELgRHARu+O7V5a/qpndvUn4k3+9SqrxD+YcUyg/p6fOfLhodGtFuhKYumarF2xZr1a5VkqT7Rtyng6GDmrJwiuZunqs4i9OU86YcM2xLUpzFafY1s7V+z3rtOrhLLy59Ud9mfquBaQM16YxJuundmzTr+1lqU7eNRrYfedhrzUwXdbtIiXGJuvDVCzXshWGKszj9fvDvi47p0qiLujTqIkn6YO0HemHpC3LOVWrEEf77aN1HMpmGth1a1Da2S0HgfuvbtzRx0ERJ0neZ32lg2kDFx8XrmcXP6I+n/1FxFqeNezdKkoa0HqJHvqzYCHeDGg1UI+Hw+eAt67RUXn6elmxborQ6aUqITtu4ddCtemDkAzr1n6fqjx//UZL017P/qpsH3Fyucy3aukjhSFg39LtBF3e7WD2a9FCcxemGjBv08NyHdeM7Nyo5IVnf7/7+sNdlh7J152l3HhaiL+t5mQa1GqSWtVsqMT6xQt+zVPDLxhfXfaFHv3xUueFcSVKjmo30h1P/oH8v/bd+8e4vdGmPS/XihS8qLz9Pd358pz5Y+4E6NOige4ffW+J7npF+huIsTh+u+7DEwJ0bztWirYt0y8BbVCOhhiZ9NkmLty7WgbwDGtdrnK496dpy1d60VlNJBdNi2jdoX+HvHThREbjhq3AkrNs/ul31Uurp7A5n+13OCeHM9meqe5PuRdvD0odJkmaum6kVO1eoYY2G6tu8r+IsTqe1Oa1S50iKT1Knhp3UqWEnndLqlMP2/X7w7zX6ldG67qTrSl2hYUyXMfrmxm+0fs96NU9tXvSx95HaN2ivfYf2KTMnU41qNqpUrfDXR99/pL7N+xaNtEpSev109WzSU+98944mDpoo55y+y/pOl/W4TINbDdaVb16pLzZ9oSGth2jj3o0ymWaNm6W9h/aW65yh/JA+XPeh5mycU3T/gFSwys66Peu0N3evXrrwJV3W87Kj6x33kb7e9rUm/GeCHv3yUd2YcWO5Rrm/3PylJGlg2kA1TW1a1J6SkKLJwyZr3LRxOrnlyfp5v58f9stjjYQa+s0pvznq/drWa1uu77U07Ru011/P+etR7eP7jVd6/XQNSy/4ZTclIUUPjHxAD4x8oMz3q5dST/1b9NfMdTN1zxn3HLV/0dZFCkVCGpQ2SPFx8Yq4iKYunipJ6teiX7nrLuy77dnbj3EkUL0QuOGrl795Wet2r9O0S6ZpTBceGloZTVObqlfTXvrnkn9q18FdOrP9mZ4uVXZep/P0xk/eOGwqSUm6Ne6mbo27lXlM4fSBtVlrCdwnoEPhQ5q3ZZ5+OeCXR+0blj5MUxZOKZqGtCd3jzo26KjRnUcrJSFF/7vsfwsC976NalG7hRLjEyt0DYzrPU7jeo+rcM0pCSk6Oe1k/eaU3+iS1y7RjNUzyvVvz5dbvlR6vfTDwnahq3pfpVEdRqlxLf+f+xAfF69RHUZV6rUj243UvXPu1d7cvap7xM2TczfNlSQNajWo6Jecf3/zb6UkpBzz57y4whHu4jeRwh8HQwf1zup3FIqUfHPukTJaZKhTw04eV1V9Ebjhiey87KJ5vmVZm7VWPZv01Pmdz49RZdXTuF7jdOfsOxVv8bqi5xWensvMdEHXC47Le7WvX/CR8rrd6wKxbnNZwpGwtu7fWqHX1Eupp9rJtT2q6NiWbFuivPy8o+YES9LQtkP12LzH9NWWr4p+AezYsKNqJ9fWeZ3O06srXtUjox7Rpr2b1Lpu61iXrgu7XqjWdVvr2reuVbOPjr1ixrrd63Rh1wtL3V8VwvaPNaLdCE36bJLu+eQedW3cVVLBVJWWtVvqw3Ufqm29tkVzsJulNtO2A9t0csuTi6btlEfBsqB21Ah3TihHn2/6/LBPLI6U0SLjsE9S8OM8s+gZ/eq9X5X7+Pb122v1hNWsPV9JBG544tMNn+q9Ne9pSOshqlPGGrXNUpvpN4N+ww/wj3TrKbfq1lNu9buMCisa4d691udK/HfFG1fo1eWvVug1DWs01JZfb1FyQvKxD/bA/C3zJanEX5ZOa3OaTKbZ62erVd2C55V1bNBRknRJ90v02orXNGfjHG3cu1EnNT8pdkVHJcQl6MlzntTzXz9fruN7Nu2pCQMmeFyVvwa1GqQmtZro4S8fLnH/Vb2uKvo6o0WGZqyeoX7Nyz+dRCro90Y1G2nbgW3aun+rvtnxjc5sf6bunn23Hvii7GkvF3e7WK9eXLGfEZRu8bbFalKriT679rNjHvv+mvf1y/d+qY+//7jEm4NxbARueGLOxjmKt3i9d8V7qlXJNW9R/dVIrKEWtVsEPnDnhHL09qq3dU7Hc3Rhl9JHUYtbsXOFHv7yYS34YYEa1Gige+fcq3xXuRuPU+JTdO+Ieyv8hNF5W+apRe0WSquTdtS+BjUaqGfTnpq9YbYGthyoOItTev10SSpa0eSzDZ9p496NGtPZn+lk53Y6V+d2OteXc1dFSfFJWvvLtdqTu0dSwZz47dnbtXX/Vjm5w5YU7de8X0HgrsD87UJNU5tqe/Z23TfnPj0+/3HN/+l8TV0yVaM6jNIfT/tjia+5//P7NXv9bG6wPo6Wbl+qPs36lGuaSOu6rXXX7Lv09KKnCdyVROCGJz7f9LlOan4SYRvH1K5+O63NCnbg/nj9x8oJ5+iXA355zLnxhXYd3KWHv3xYn274VOt2r9Mry15Rm3ptKnzu/Ei+vt/zvU5pdYqu73t9hV47b8u8o9aJLm5om6F6etHTCkfCalO3TdFScfVS6qlb4256e/XbOpR/yJcpJShZalLqYetjF346caTh6cP158/+rFNbl7yEYFma1mqq7Qe2a3fObknS6FdGa9fBXZo4cOJRN2kXOr/T+Zq+arpWZ64u9SZslF84EtbyHcvL/alNSkKKrup1lf6+8O+a9OmkSn0qHWdxuqrXVWpZp2WFX1sdELhx3OXl52nelnm6od8NfpeCE0D7+u314boP/S7DVzNWz1CtxFo6ve3p5X5No5qN1K1xN83eMFtLti3Rf3X7L738Xy9X+NwRF1Hqn1MrvB56Vk6W1mSt0fUnlR7Sz+pwlh6f/7hmr5+tsV3GHrZvYMuBmrqkYJULAveJ59Q2p2rvbXtVI7HGsQ8+QtPUppq7aa72Htqrmok1te3ANqXXSz9sLfcjDWk9RFLBp6cEbmlv7l79Ze5fij6NKE28xeu6k65Tz6Y9D2v/LvM7Hco/pF5Ne5X7nDdk3KCnFz1dtKxmZby16i19cd0XgfyUgsCN427x1sXKDedqcOvSn2gIFGpfv72e3/+8ckI5lfrP+0TnnNOM1TM0sv3IMh8PXpJTW5+qKQunyMlpbOexlTp/nMWpW+NuWr5zeYVeVzR/u2XpN7ue3eFsLbtxmXLDuUeFpEGtBhUF7tJGUVG1VfbntVmtZlq/Z72cnCadMUnPLn5Wtw66tcxR084NO6thjYaas2lOhT+JOd7CkbB2ZO/w7fwHQwd12euXadHWRaqbXPeYxz67+Fm99pPXdGb7M4val25fKkkVCtxdG3fV/tv3V3rq2nNLntPPZ/xcb377Zpk3IFdXBG4cZtWuVZqycIrOSD9Dreu21qHwIeXl55V55/iRCh/zW9YjxIFChQ/jqX9//UCOejjndCj/kO48/c4Kv/a0NqfpHwv/ocS4RJ3dsfLr2Hdv0l0frq3YpwxfbPpCcRZX5hxeMztszfjiiq9swgh3sDRNbSonJ0ka3Hqw7jjtjmO+xsw0pPUQzdk4x+vyyrQ2a63Ofelcrcpc5WsdSfFJmnbJtGOu8LV532ad+9K5Gv3yaH19w9dFv/gu3b5UCXEJRQ8hK6/4uHjFq3xPZz3SdSddp0e/fFQT/jNBLy+r+KdxR4qzOHVq0EkdG3Y86pe13k17HzWq7zcCNw5z6we36p3v3in1LvXy6tyws5rXbn6cqkJ1dl6n83TX6XcpJ5Tjdym+qZFYQ5f1OPohLcdSOH92WPqwMlcDOpYejXvoha9fUFZOVrmXXftw3Yca0HJApc/brXE31Umuo1B+SA1rNKzUe+DEVLgWtyT1bFL+UHRq61P11qq3NPCZgb79cr5q1yqZmR4961FfP5Eb0HKA+jTrc8zj0uqk6f0r31fXJ7rqp2//VJ9c84niLE5LdyxVl0ZdYrrCUUJcgqacP0W/eu9XWrFzxY9+v1B+SK+veL3EEfd7ht5D4EbV9fW2r/XOd+/oj6f9Uae3OV17cvcoOSFZSfFJireK/UbL4vgor9rJtXX30Lv9LuOE1KpuK/1m0G9+9EobPZr0kCQt37G8xMd6H2lP7h7N3zJfd5x67JHJ0sRZnAalDdLmfZsD+clGkBU+PKh5anM1rFn+X7Yu7n6xPt34adHj7P0wLH2Y/jz8zyfU/3HNUpvpL2f+RddPv15tHm2jxLhEbd63WRd3vzjmtQxpPUQLxy88bu93MHRQP+z/4aj2qrheO4Eb2pG9Q9NXTdcry15RalKqJg6cqPo16vtdFoByePDMB3/0exRO+1i+s3yB++PvP1bERTSy3cgfdd5/nPcPHcg78KPeAyeewhHuio5Atq7bWm9d+pYXJVV71/a5VjuydxSNLJuZfpHxC5+r+vFqJtYsmpZY1RG4A27Lvi0a+vxQrclaI0m649Q7CNtAwLSq00q1k2pr+Y7y3Tj54boPlZqUWuITJiuiMssY4sRX+LTKikwnwY9jZrptyG1+lxFoBO4Aygnl6M6P79TrK19XZk6mnHOaedVMdW/S/bC5dQCCofDmxtdXvn7UI7dLMuv7WRradqgS4xNjUB2qm2apzXTroFt1TZ9r/C4FiBkCdzX39wV/12PzHjusLSsnSzuyd+j8TuerUc1GujHjRvVv2d+nCgFUBeN6jdNf5/+1XOtxN0ttpvF9x8egKlRHZqaHznzI7zKAmDLnnN81eCojI8MtWLDA7zJ80/3J7jqQd+Cwj34T4hJ0bZ9ry3zIAAAAAMrPzBY65zJK2scIdzX2w/4ftGLnCj0w4gH9dvBv/S4HAAAgkEp/rBNOeLO+nyVJjGQDAAD4iBHuE8xtM2/TzHUzi7av7n21Jpw8ocRjZ66bqYY1Gqp3s96xKg8AAABHYIT7BOKc09/m/017D+1Vs9Rm2pO7R3fNvqvEhwA45/TR9x9pWPqwox55CgAAgNghiZ1ANu3bpOxQtm4ddKtmXD5DT577pHbn7tb0VdOPOnbu5rnavG+zhqcP96FSAAAAFCJwn0BW7lwpSeraqKskaXj6cKXVSdM/l/zzsONyw7n66fSfKq1Omi7reVnM6wQAAMD/YQ53FbE3d6+G/HOIVuxcofb122v5L5Yf9VCJlbuigbtxQeCOj4vX1b2v1r1z7tWQqUOKjtudu1srd63Ue1e8pzrJdWL3TQAAAOAojHBXEZM/m6zlO5ZrbJex+i7rOy3ZtuSoY1buXKmGNRqqcc3GRW03Ztyoczqeo5SElKI/zVOb66GRD+msDmfF8DsAAABASRjhrgLW7V6nx+Y9pqv7XK3JwybrjZVvaM7GOUc9/XHlrpXq2rirzKyorWWdlnr7srdjXTIAAADKiRHuKuCv8/4qk2nysMlqUbuF2tVvpzmb5hx13IqdK4rmbwMAAODEQOCuAj7Z8IlOaXWKWtRuIUka0nqI5mycI+dc0TE7s3cqMyeTwA0AAHCCIXD7bN+hffp6+9c6tfWpRW1DWg3RjuwdWpO1pqit8IbJbo27xbxGAAAAVB5zuH32xaYvFHERndqmWOBuXbDiyF2z71KPJj2UE8rRM4ufUVJ8Ek+NBAAAOMEQuH322YbPFG/xGpg2sKitS6Mu6tKoi15e9nJR25DWQ/TgyAfVLLWZH2UCAACgkgjcPvt046fq27yvUpNSi9rMTMt/sVzhSLioLSk+yY/yAAAA8CMxh9snL3z9glo90kpfbPrisPnbheIsTknxSUV/AAAAcGIicPtk9vrZ2p2zWz/r+zPd2P9Gv8sBAACAR6pc4DazB83sWzNbamZvmlm9YvtuN7M1ZrbKzE7oxyhmh7KVVidNfz/v7+rQoIPf5QAAAMAjVS5wS/pQUg/nXC9JqyXdLklm1k3SpZK6Sxol6Ukzi/etyh8pOy9btZJq+V0GAAAAPFblArdz7gPnXOHdgl9KSot+PUbSK865Q8657yWtkTTAjxqPh+xQtmolErgBAACquyoXuI9wnaT/RL9uKWlTsX2bo21HMbPxZrbAzBbs3LnT4xIrhxFuAACAYPAlcJvZTDNbVsKfMcWOuUNSWNKLFX1/59wU51yGcy6jcePGx7P044YRbgAAgGDwZR1u59yIsvab2TWSzpM03Dnnos1bJLUqdlhatO2EdCDvACPcAAAAAVDlppSY2ShJv5M02jl3sNiu6ZIuNbNkM0uX1FHSfD9qPB6y8xjhBgAACIKq+KTJv0lKlvShmUnSl865G5xzy83sVUkrVDDV5CbnXL6Pdf4oTCkBAAAIhioXuJ1zpS5K7ZybLGlyDMvxRH4kX7nhXKaUAAAABECVm1ISBAdDBTNlGOEGAACo/gjcPsgOZUsSI9wAAAABQOD2QXZeNHAzwg0AAFDtEbh9wAg3AABAcBC4fcAINwAAQHAQuH3ACDcAAEBwELh9wAg3AABAcBC4fcAINwAAQHAQuH3ACDcAAEBwELh9wAg3AABAcBC4fcAINwAAQHAQuH2QHcpWvMUrKT7J71IAAADgMQK3D7LzslUrqZbMzO9SAAAA4DECtw+yQ9lMJwEAAAgIArcPskPZ3DAJAAAQEARuH2TnMcINAAAQFARuHzDCDQAAEBwEbh8wwg0AABAcBG4fMMINAAAQHARuHzDCDQAAEBwEbh+wLCAAAEBwELh9UPjgGwAAAFR/BO4Yc87pYOggI9wAAAABQeCOsZxwjpwcI9wAAAABkeB3AUHyzfZv9K+l/5IkRrgBAAACgsAdQw9+8aD+tfRfSo5PVvcm3f0uBwAAADFA4I6hQ/mH1KVRF628aaXfpQAAACBGmMMdQ6H8kBLi+B0HAAAgSAjcMRSOhAncAAAAAUPgjqFQJKTEuES/ywAAAEAMEbhjiBFuAACA4CFwx1AoP6TEeEa4AQAAgoTAHUOMcAMAAAQPgTuGmMMNAAAQPATuGGKEGwAAIHiqbOA2s1vNzJlZo+i2mdnjZrbGzJaaWV+/a6yocCTMHG4AAICAqZKB28xaSTpT0sZizWdL6hj9M17SUz6U9qPw4BsAAIDgqZKBW9Ijkn4nyRVrGyPpBVfgS0n1zKy5L9VVUjgSZg43AABAwFS5wG1mYyRtcc59fcSulpI2FdveHG0r6T3Gm9kCM1uwc+dOjyqtuFCEEW4AAICg8SX9mdlMSc1K2HWHpD+oYDpJpTnnpkiaIkkZGRnuGIfHDDdNAgAABI8v6c85N6KkdjPrKSld0tdmJklpkhaZ2QBJWyS1KnZ4WrTthBHKZ1lAAACAoKlSU0qcc98455o459o659qqYNpIX+fcNknTJY2LrlYyUNJe59xWP+utKEa4AQAAgudESn/vSjpH0hpJByVd6285FReK8Gh3AACAoKnSgTs6yl34tZN0k3/V/HiMcAMAAARPlZpSUt2xLCAAAEDwELhjiAffAAAABA+BO0YiLiInxxxuAACAgCFwx0goPyRJjHADAAAEDIE7RsKRsCQCNwAAQNAQuGMkFCkY4eamSQAAgGAhcMcII9wAAADBROCOkcI53Nw0CQAAECwE7hhhhBsAACCYCNwxwhxuAACAYCJwxwgj3AAAAMFE4I6RwsDNHG4AAIBgIXDHCA++AQAACCYCd4wUjXAzhxsAACBQCNwxUnjTJCPcAAAAwULgjhFumgQAAAgmAneM8OAbAACAYCJwxwgj3AAAAMFE4I4RHnwDAAAQTATuGGGEGwAAIJgI3DHCg28AAACCicAdIzz4BgAAIJgI3DHCg28AAACCicAdIzz4BgAAIJgI3DHCHG4AAIBgInDHCHO4AQAAgqnU9Gdmb0type13zo32pKJqimUBAQAAgqms9PdQ9O8LJTWT9O/o9mWStntZVHXEg28AAACCqdTA7Zz7RJLM7C/OuYxiu942swWeV1bNMMINAAAQTOWZw13LzNoVbphZuqRa3pVUPRXO4eamSQAAgGApz3DrREmzzWydJJPURtJ4T6uqhhjhBgAACKZjpj/n3Htm1lFSl2jTt865Q96WVf2EI2HFWZzijIVhAAAAguSYgdvMEiX9XNJp0abZZvYP51zI08qqmVAkxOg2AABAAJUnAT4lKVHSk9Htq6JtP/WqqOooHAmzQgkAAEAAlSdw93fO9S62PcvMvvaqoOoqlM8INwAAQBCVZ0Jxvpm1L9yIrliS711JkplNMLNvzWy5mT1QrP12M1tjZqvM7CwvazjewpEwK5QAAAAEUHmGXH8r6eMjVim51quCzOwMSWMk9XbOHTKzJtH2bpIuldRdUgtJM82sk3PO0/B/vDCHGwAAIJjKs0rJR9FVSjpHm1Z5vErJjZLuKzyHc25HtH2MpFei7d+b2RpJAyTN9bCW4yYcCRO4AQAAAuiYU0qKrVJyZ/TPz6JtXukk6VQzm2dmn5hZ/2h7S0mbih23OdpWUs3jzWyBmS3YuXOnh6WWXygS4qZJAACAAPJllRIzmympWQm77ojW1EDSQEn9Jb1a/EmX5eGcmyJpiiRlZGS4ytZ5PDHCDQAAEEy+rFLinBtR2j4zu1HSG845J2m+mUUkNZK0RVKrYoemRdtOCNw0CQAAEExVcZWSaZLOiJ6rk6QkSbskTZd0qZklm1m6pI6S5ntYx3HFsoAAAADBVOVWKZE0VdJUM1smKU/S1dHR7uVm9qqkFZLCkm46UVYokXjwDQAAQFBVuVVKnHN5kq4sZd9kSZO9OreXWBYQAAAgmMqbAPtJahs9vo+ZyTn3gmdVVUPM4QYAAAimYwZuM/uXpPaSluj/5m47SQTuCmAONwAAQDCVJwFmSOoWnUeNSgpHwkpJSPG7DAAAAMRYeVYpWaaS18xGBTCHGwAAIJhKTYBm9rYKpo7UlrTCzOZLKrpZ0jk32vvyqg8efAMAABBMZSXAh2JWRQCE8kPcNAkAABBApQZu59wnsSykumOEGwAAIJjKmlIyxzk3xMz2q2BqSdEuSc45V8fz6qoRHnwDAAAQTGWNcA+J/l07duVUX9w0CQAAEExljXA3KOuFzrms419O9cUINwAAQDCVNeS6UAVTSayEfU5SO08qqqZ48A0AAEAwlTWlJD2WhVR3PNodAAAgmI754BsrcKWZ/TG63drMBnhfWvXCHG4AAIBgKs+TJp+UNEjS5dHt/ZKe8Kyiaoo53AAAAMFUniHXk51zfc1ssSQ553abWZLHdVU7zOEGAAAIpvKMcIfMLF7RtbjNrLGkiKdVVTPOOeW7fAI3AABAAJUncD8u6U1JTcxssqQ5kv7saVXVTL7LlyRumgQAAAig8gy5vqaCJQKHq2CJwLGStntYU7UTyg9JEiPcAAAAAVSeBPiGpLHOuW8lycyaS/pQUj8vC6tOwpGwJHHTJAAAQACVZ0rJNEmvmlm8mbWV9L6k270sqroJRRjhBgAACKpjJkDn3NPRVUmmSWor6efOuS88rqtaKRrhZg43AABA4JQauM3s18U3JbWWtETSQDMb6Jx72OPaqg3mcAMAAARXWQmw9hHbb5TSjmNgDjcAAEBwlRq4nXN/imUh1RlzuAEAAIKrrCkljzrnbjGztxV96E1xzrnRnlZWjRSOcBO4AQAAgqesBPiv6N8PxaKQ6oybJgEAAIKrrCklC6N/fxK7cqonRrgBAACCq6wpJd+ohKkkhZxzvTypqBoicAMAAARXWQnwvJhVUc2xLCAAAEBwlTWlZEMsC6nOGOEGAAAIrvI82h0/EutwAwAABBeBOwYY4QYAAAguAncMELgBAACC65gJsJTVSvZKWiBpknMu04vCqhMCNwAAQHCVJwH+R1K+pJei25dKqilpm6TnJJ3vSWXVCIEbAAAguMqTAEc45/oW2/7GzBY55/qa2ZXHuyAz6yPp75JSJIUl/cI5N9/MTNJjks6RdFDSNc65Rcf7/F4gcAMAAARXeeZwx5vZgMINM+svKT66Gfagpgck/ck510fSndFtSTpbUsfon/GSnvLg3J4IRViHGwAAIKjKkwB/KmmqmaVKMkn7JF1vZrUk3etBTU5SnejXdSX9EP16jKQXnHNO0pdmVs/MmjvntnpQw3HFCDcAAEBwHTMBOue+ktTTzOpGt/cW2/2qBzXdIul9M3tIBSPwp0TbW0raVOy4zdG2owK3mY1XwSi4Wrdu7UGJFVO0Dnc863ADAAAETXlWKakr6S5Jp0W3P5F0zxHBu0LMbKakZiXsukPScEkTnXOvm9lPJD0raURF3t85N0XSFEnKyMg4coWVmGOEGwAAILjKkwCnSlom6SfR7ask/VPShZU9qXOu1ABtZi9I+lV08/9Jeib69RZJrYodmhZtq/II3AAAAMFVnpsm2zvn7nLOrYv++ZOkdh7W9IOk06NfD5P0XfTr6ZLGWYGBkvaeCPO3JQI3AABAkJUnAeaY2RDn3BxJMrPBknI8rOlnkh4zswRJuYrOxZb0rgqWBFyjgmUBr/WwhuOKwA0AABBc5UmAN0h6ofCmSUm7JV3tVUHRYN+vhHYn6SavzuulUD7LAgIAAARVeVYp+VpSbzOrE93eZ2a3SFrqcW3VBiPcAAAAwVWeOdySCoK2c25fdPPXHtVTLRG4AQAAgqvcgfsIdlyrqObCkbDiLE5xVtnuBgAAwImqsgnQ97WtTyThSJjRbQAAgIAqNQWa2X6VHKxNUg3PKqqGCNwAAADBVWoKdM7VjmUh1RmBGwAAILiYVBwDBG4AAIDgInDHQCgSInADAAAEFIE7BhjhBgAACC4CdwyEI2ElxiX6XQYAAAB8QOCOAUa4AQAAgovAHQMEbgAAgOAicMcAgRsAACC4CNwxQOAGAAAILgJ3DLAsIAAAQHARuGOAEW4AAIDgInDHQDgSVmI8ywICAAAEEYE7BhjhBgAACC4CdwwQuAEAAIKLwB0DBG4AAIDgInDHAIEbAAAguAjcMUDgBgAACC4CdwyE8lmHGwAAIKgI3DHACDcAAEBwEbhjIBwJKzGOdbgBAACCiMAdA4xwAwAABBeBOwYI3AAAAMFF4I4BAjcAAEBwEbhjgMANAAAQXATuGAhFWBYQAAAgqAjcMcAINwAAQHARuGOAwA0AABBcBG6POedYhxsAACDACNwei7iIJDHCDQAAEFC+BG4zu9jMlptZxMwyjth3u5mtMbNVZnZWsfZR0bY1ZnZb7KuunHAkLInADQAAEFR+jXAvk3ShpE+LN5pZN0mXSuouaZSkJ80s3sziJT0h6WxJ3SRdFj22yiNwAwAABJsvKdA5t1KSzOzIXWMkveKcOyTpezNbI2lAdN8a59y66OteiR67IjYVVx6BGwAAINiq2hzulpI2FdveHG0rrb3KC0VCkgjcAAAAQeVZCjSzmZKalbDrDufcW16dN3ru8ZLGS1Lr1q29PNUxMcINAAAQbJ6lQOfciEq8bIukVsW206JtKqO9pHNPkTRFkjIyMlwl6jhuCgN3YjzLAgIAAARRVZtSMl3SpWaWbGbpkjpKmi/pK0kdzSzdzJJUcGPldB/rLDdGuAEAAILNlxRoZhdI+qukxpLeMbMlzrmznHPLzexVFdwMGZZ0k3MuP/qamyW9Lyle0lTn3HI/aq8oAjcAAECw+bVKyZuS3ixl32RJk0tof1fSux6XdtwRuAEAAIKtqk0pqXYI3AAAAMFG4PZYKJ9lAQEAAIKMwO0xRrgBAACCjcDtMQI3AABAsBG4PVa0Dncc63ADAAAEEYHbY4xwAwAABBuB22MEbgAAgGAjcHuMwA0AABBsBG6PEbgBAACCjcDtsVCEdbgBAACCjMDtMUa4AQAAgo3A7bGiZQHjWRYQAAAgiAjcHmOEGwAAINgI3B4jcAMAAAQbgdtjBG4AAIBgI3B7jMANAAAQbARuj4XyWRYQAAAgyAjcHmOEGwAAINgI3B4jcAMAAAQbgdtjRetwx7EONwAAQBARuD1WGLjj4+J9rgQAAAB+IHB7LBwJK87iFGd0NQAAQBCRAj0WjoSZvw0AABBgBG6PEbgBAACCjcDtsVAkROAGAAAIMAK3xxjhBgAACDYCt8cI3AAAAMFG4PZYOBJmDW4AAIAAI3B7jBFuAACAYCNwe4zADQAAEGwEbo8RuAEAAIKNwO0xlgUEAAAINgK3x/Ly85QYz02TAAAAQUXg9tjunN2qn1Lf7zIAAADgEwK3x7JystSgRgO/ywAAAIBPfAncZnaxmS03s4iZZRRrH2lmC83sm+jfw4rt6xdtX2Nmj5uZ+VF7RWXlZKlhjYZ+lwEAAACf+DXCvUzShZI+PaJ9l6TznXM9JV0t6V/F9j0l6WeSOkb/jIpBnT+Kc06ZOZmMcAMAAASYL8tnOOdWStKRg9TOucXFNpdLqmFmyZIaSKrjnPsy+roXJI2V9J9Y1FtZB/IOKBwJq2FNRrgBAACCqirP4f4vSYucc4cktZS0udi+zdG2Ki0rJ0uSGOEGAAAIMM9GuM1spqRmJey6wzn31jFe213S/ZLOrOS5x0saL0mtW7euzFscF5k5mZII3AAAAEHmWeB2zo2ozOvMLE3Sm5LGOefWRpu3SEordlhatK20c0+RNEWSMjIyXGXqOB4KR7i5aRIAACC4qtSUEjOrJ+kdSbc55z4vbHfObZW0z8wGRlcnGSepzFHyqiDzICPcAAAAQefXsoAXmNlmSYMkvWNm70d33Sypg6Q7zWxJ9E+T6L5fSHpG0hpJa1XFb5iUmMMNAAAA/1YpeVMF00aObJ8kaVIpr1kgqYfHpR1XBG4AAABUqSkl1U1mTqZqJdZSckKy36UAAADAJwRuD2XlZLEGNwAAQMARuD3EUyYBAABA4PZQVk4WgRsAACDgCNweysrJYg1uAACAgCNweyjzIFNKAAAAgo7A7RHnHCPcAAAAIHB7ZX/efuW7fEa4AQAAAs6XB98EAY91BwAAXguFQtq8ebNyc3P9LiUwUlJSlJaWpsTExHK/hsDtkcKnTLIONwAA8MrmzZtVu3ZttW3bVmbmdznVnnNOmZmZ2rx5s9LT08v9OqaUeCQzhxFuAADgrdzcXDVs2JCwHSNmpoYNG1b4EwUCt0cKR7gJ3AAAwEuE7diqTH8TuD1SNKWEVUoAAEA1N3nyZHXv3l29evVSnz59NG/ePM/ONXToUC1YsMCz9/cCc7g9UnjTZP0a9X2uBAAAwDtz587VjBkztGjRIiUnJ2vXrl3Ky8vzu6wqhRFuj2TlZKl2Um0lxSf5XQoAAIBntm7dqkaNGik5OVmS1KhRI7Vo0UL33HOP+vfvrx49emj8+PFyzkkqGKGeOHGiMjIy1LVrV3311Ve68MIL1bFjR/33f/+3JGn9+vXq0qWLrrjiCnXt2lUXXXSRDh48eNS5P/jgAw0aNEh9+/bVxRdfrAMHDsTuG68ARrg9kpWbxfxtAAAQM7e8d4uWbFtyXN+zT7M+enTUo2Uec+aZZ+qee+5Rp06dNGLECF1yySU6/fTTdfPNN+vOO++UJF111VWaMWOGzj//fElSUlKSFixYoMcee0xjxozRwoUL1aBBA7Vv314TJ06UJK1atUrPPvusBg8erOuuu05PPvmkfvOb3xSdd9euXZo0aZJmzpypWrVq6f7779fDDz9cdM6qhBFuj/BYdwAAEASpqalauHChpkyZosaNG+uSSy7Rc889p48//lgnn3yyevbsqVmzZmn58uVFrxk9erQkqWfPnurevbuaN2+u5ORktWvXTps2bZIktWrVSoMHD5YkXXnllZozZ85h5/3yyy+1YsUKDR48WH369NHzzz+vDRs2xOi7rhhGuD2SlZPFGtwAACBmjjUS7aX4+HgNHTpUQ4cOVc+ePfWPf/xDS5cu1YIFC9SqVSvdfffdhy2lVzj9JC4urujrwu1wOCzp6NVAjtx2zmnkyJF6+eWXvfq2jhtGuD2SmcMINwAAqP5WrVql7777rmh7yZIl6ty5s6SC+dwHDhzQa6+9VuH33bhxo+bOnStJeumllzRkyJDD9g8cOFCff/651qxZI0nKzs7W6tWrK/tteIoRbo9k5WSxJCAAAKj2Dhw4oAkTJmjPnj1KSEhQhw4dNGXKFNWrV089evRQs2bN1L9//wq/b+fOnfXEE0/ouuuuU7du3XTjjTcetr9x48Z67rnndNlll+nQoUOSpEmTJqlTp07H5fs6nqzwjtHqKiMjw8V6rcaIiyjxfxJ1+5DbNWnYpJieGwAABMfKlSvVtWtXv8s47tavX6/zzjtPy5Yt87uUEpXU72a20DmXUdLxTCnxwL5D+xRxEaaUAAAAgMDtBZ4yCQAAUHlt27atsqPblUHg9kDhUyYZ4QYAAACB2wNFI9wsCwgAABB4BG4PFAZuRrgBAABA4PZAZg5TSgAAAFCAwO0BRrgBAEBQbN++XZdffrnatWunfv36adCgQXrzzTc9Odfdd9+thx566Kj2adOmacWKFRV+v/Xr1+ull14q2n7uued08803/6gaS0Lg9kBWTpbqJNdRQhzPFQIAANWXc05jx47VaaedpnXr1mnhwoV65ZVXtHnz5qOOLXxkuxfKCtxlnffIwO0VArcHeKw7AAAIglmzZikpKUk33HBDUVubNm00YcIESQUjxqNHj9awYcM0fPhwZWVlaezYserVq5cGDhyopUuXSioYub7uuus0dOhQtWvXTo8//njR+02ePFmdOnXSkCFDtGrVqqNq+OKLLzR9+nT99re/VZ8+fbR27VoNHTpUt9xyizIyMvTYY4/pmmuuOezx8qmpqZKk2267TZ999pn69OmjRx55RJL0ww8/aNSoUerYsaN+97vfHZd+YgjWAzzWHQAAxNwtt0hLlhzf9+zTR3r00VJ3L1++XH379i3zLRYtWqSlS5eqQYMGmjBhgk466SRNmzZNs2bN0rhx47QkWvO3336rjz/+WPv371fnzp114403aunSpXrllVe0ZMkShcNh9e3bV/369Tvs/U855RSNHj1a5513ni666KKi9ry8PBU+bfyaa64psbb77rtPDz30kGbMmCGp4BeEJUuWaPHixUpOTlbnzp01YcIEtWrVqux+OgYCtwd6NO6hzg07+10GAABATN10002aM2eOkpKS9NVXX0mSRo4cqQYNCj75nzNnjl5//XVJ0rBhw5SZmal9+/ZJks4991wlJycrOTlZTZo00fbt2/XZZ5/pggsuUM2aNSVJo0ePLnctl1xySaW+h+HDh6tu3bqSpG7dumnDhg0E7qro/pH3+10CAAAImjJGor3SvXv3ogAtSU888YR27dqljIyMorZatWqV672Sk5OLvo6Pj//Rc76LnzchIUGRSESSFIlElJeXF7M6JOZwAwAAoJKGDRum3NxcPfXUU0VtBw8eLPX4U089VS+++KIkafbs2WrUqJHq1KlT6vGnnXaapk2bppycHO3fv19vv/12icfVrl1b+/fvL/V92rZtq4ULF0qSpk+frlAoVK7XHS8EbgAAAFSKmWnatGn65JNPlJ6ergEDBujqq6/W/feX/Gn/3XffrYULF6pXr1667bbb9Pzzz5f5/n379tUll1yi3r176+yzz1b//v1LPO7SSy/Vgw8+qJNOOklr1649av/PfvYzffLJJ+rdu7fmzp1bNPrdq1cvxcfHq3fv3kU3TXrBnHOevXmpJzW7WNLdkrpKGuCcW3DE/taSVki62zn3ULRtlKTHJMVLesY5d195zpWRkeEKJ8wDAABUJytXrlTXrl39LiNwSup3M1vonMso6Xi/RriXSbpQ0qel7H9Y0n8KN8wsXtITks6W1E3SZWbWzesiAQAAgB/Ll5smnXMrpYKPIY5kZmMlfS8pu1jzAElrnHProse8ImmMCkbBAQAAgCqrSs3hNrNUSb+X9KcjdrWUtKnY9uZoGwAAAFCleTbCbWYzJTUrYdcdzrm3SnnZ3ZIecc4dKGn0uwLnHi9pvCS1bt260u8DAABQ1TnnSpw1AG9U5v5HzwK3c25EJV52sqSLzOwBSfUkRcwsV9JCScVXHE+TtKWMc0+RNEUquGmyEnUAAABUeSkpKcrMzFTDhg0J3THgnFNmZqZSUlIq9Loq9eAb59yphV+b2d2SDjjn/mZmCZI6mlm6CoL2pZIu96dKAACAqiEtLU2bN2/Wzp07/S4lMFJSUpSWllah1/gSuM3sAkl/ldRY0jtmtsQ5d1ZpxzvnwmZ2s6T3VbAs4FTn3PLYVAsAAFA1JSYmKj093e8ycAy+rMMdS6zDDQAAAK9VxXW4AQAAgEAgcAMAAAAeqvZTSsxsp6QNPpy6kaRdPpz3REV/VRx9VjH0V8XQXxVDf1UcfVYx9FfF+NFfbZxzjUvaUe0Dt1/MbEFp83hwNPqr4uiziqG/Kob+qhj6q+Los4qhvyqmqvUXU0oAAAAADxG4AQAAAA8RuL0zxe8CTjD0V8XRZxVDf1UM/VUx9FfF0WcVQ39VTJXqL+ZwAwAAAB5ihBsAAADwEIHbA2Y2ysxWmdkaM7vN73qqIjNbb2bfmNkSM1sQbWtgZh+a2XfRv+v7XadfzGyqme0ws2XF2krsHyvwePR6W2pmff2r3B+l9NfdZrYleo0tMbNziu27Pdpfq8zsLH+q9o+ZtTKzj81shZktN7NfRdu5xkpRRp9xnZXAzFLMbL6ZfR3trz9F29PNbF60X/7XzJKi7cnR7TXR/W19/QZirIz+es7Mvi92ffWJtgf+Z1KSzCzezBab2YzodpW9vgjcx5mZxUt6QtLZkrpJuszMuvlbVZV1hnOuT7Fle26T9JFzrqOkj6LbQfWcpFFHtJXWP2dL6hj9M17SUzGqsSp5Tkf3lyQ9Er3G+jjn3pWk6M/jpZK6R1/zZPTnNkjCkm51znWTNFDSTdF+4RorXWl9JnGdleSQpGHOud6S+kgaZWYDJd2vgv7qIGm3pOujx18vaXe0/ZHocUFSWn9J0m+LXV9Lom38TBb4laSVxbar7PVF4D7+Bkha45xb55zLk/SKpDE+13SiGCPp+ejXz0sa618p/nLOfSop64jm0vpnjKQXXIEvJdUzs+YxKbSKKKW/SjNG0ivOuUPOue8lrVHBz21gOOe2OucWRb/er4L/sFqKa6xUZfRZaQJ9nUWvlQPRzcToHydpmKTXou1HXmOF195rkoabmcWmWv+V0V+lCfzPpJmlSTpX0jPRbVMVvr4I3MdfS0mbim1vVtn/KAeVk/SBmS00s/HRtqbOua3Rr7dJaupPaVVWaf3DNVe6m6Mft061/5uiRH8VE/1o9SRJ88Q1Vi5H9JnEdVai6Mf9SyTtkPShpLWS9jjnwtFDivdJUX9F9++V1DCmBfvsyP5yzhVeX5Oj19cjZpYcbQv89SXpUUm/kxSJbjdUFb6+CNzwyxDnXF8VfCx2k5mdVnynK1g+hyV0SkH/lMtTktqr4OPZrZL+4ms1VZCZpUp6XdItzrl9xfdxjZWshD7jOiuFcy7fOddHUpoKRve7+FtR1XZkf5lZD0m3q6Df+ktqIOn3/lVYdZjZeZJ2OOcW+l1LeRG4j78tkloV206LtqEY59yW6N87JL2pgn+Mtxd+JBb9e4d/FVZJpfUP11wJnHPbo/+BRSQ9rf/7OJ/+kmRmiSoIji86596INnONlaGkPuM6Ozbn3B5JH0sapIKpDwnRXcX7pKi/ovvrSsqMbaVVQ7H+GhWdyuScc4ck/VNcX4UGSxptZutVMHV3mKTHVIWvLwL38feVpI7RO2WTVHDTzHSfa6pSzKyWmdUu/FrSmZKWqaCfro4edrWkt/ypsMoqrX+mSxoXvWt9oKS9xaYFBNYR8xkvUME1JhX016XRu9bTVXDT0fxY1+en6NzFZyWtdM49XGwX11gpSuszrrOSmVljM6sX/bqGpJEqmPf+saSLoocdeY0VXnsXSZrlAvSgkFL669tivwCbCuYjF7++Avsz6Zy73TmX5pxrq4KcNcs5d4Wq8PWVcOxDUBHOubCZ3SzpfUnxkqY655b7XFZV01TSm9H7FRIkveSce8/MvpL0qpldL2mDpJ/4WKOvzOxlSUMlNTKzzZLuknSfSu6fdyWdo4Kbsg5KujbmBfuslP4aGl1Cy0laL+nnkuScW25mr0paoYKVJ25yzuX7ULafBku6StI30TmjkvQHcY2VpbQ+u4zrrETNJT0fXZklTtKrzrkZZrZC0itmNknSYhX8EqPo3/8yszUquAH6Uj+K9lFp/TXLzBpLMklLJN0QPZ6fyZL9XlX0+uJJkwAAAICHmFICAAAAeIjADQAAAHiIwA0AAAB4iMANAAAAeIjADQAAAHiIwA0AAAB4iMANAAAAeIjADQAAAHjo/wMXQGN8Sf/XOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(ppcs, label=\"Sample\", c=\"g\")\n",
    "# plotting the ground truth for reference\n",
    "plt.plot(\n",
    "    [log_likelihood(x_obs, y_obs, thetas_obs, mus_obs, sigmas_obs, N)] * num_samples,\n",
    "    label=\"Grond truth\",\n",
    "    c=\"r\",\n",
    ")\n",
    "plt.ylabel(\"Log likelihood\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "customInput": null,
    "originalKey": "6d0a7b07-fb2f-447c-af73-4e0cb86f0a84",
    "showInput": false
   },
   "source": [
    "From the above plot, inference appears to be doing a good job of fitting the random\n",
    "variables given the observed data. Inference appears to converge with a log likelihood\n",
    "scores near to those generated by the ground truth parameters."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
